| 查看: 2628 | 回复: 5 | |||
| 【悬赏金币】回答本帖问题，作者dsmj1995将赠送您 10 个金币 | |||
[求助]
请教如何提高MATLAB函数norminv(p,0,1)的使用精度
|
|||
|
请教各位大佬一个MATLAB问题：如何提高函数norminv(p,0,1)的使用精度？ norminv为一维正态分布的累积分布函数的逆函数，假设 p = 1-10^(i) , i = 1,...,.N， 因为matlab的运算精度为双精度16位有效数字 ，因此当N>16的时候，norminv(p,0,1)返回的都是正无穷大。现在编程涉及到的p比较小，所以想请教一下如何使得p<1-10^16时，也能使得该函数得出较为准确的结果？ |
回复此楼» 猜你喜欢
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有2人回复
-
070200求调剂，一志愿某211，288分
已经有18人回复
-
物理学I论文润色/翻译怎么收费?
已经有289人回复
-
法国博士后职位
已经有0人回复
-
重庆交大26年硕士生招生拟调剂通知已出！欢迎加入机器视觉与3D光学成像课题组。
已经有0人回复
-
**
已经有1人回复
-
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有0人回复
-
基底STO,薄膜SRO，XRD里面的振荡，是laue震荡还是kiessig振荡? 怎么判断？
已经有0人回复
kingspin
木虫 (著名写手)
- 应助: 254 (大学生)
- 金币: 4472.5
- 红花: 74
- 帖子: 1635
- 在线: 315.5小时
- 虫号: 2585875
- 注册: 2013-08-07
- 性别: GG
- 专业: 聚合物共混与复合材料
2楼2020-08-29 20:16:43
|
该问题已经解决。实际上，我本来的目的是为生成服从左截断的正态分布变量，X~N(mu,sigma,u-), 式中mu,sigma为标准正态变量的均值与方差，u-为左截断的位置。通常的的方法是生成均匀分布变量 p~U(p1,1)，其中p1 = normcdf(u-,mu,sigma)，然后再根据正态函数的累积分布函数将p1转化目标的截断正态分布变量X， 即是X = norminv(p1,mu,sigma)。其中normcdf，norminv都是MATLAB内置函数，分别是正态分布函数的累积分布函数与其逆函数。这种做法在u-与mu比较接近的时候，效率还是可以的。但是当u->>mu，即是p1非常接近于1，如同我在帖子中所描述的问题，由于数值计算的精度所限制这种做法是行不通的。 关于当p1非常接近于1时，对应的逆函数值的求解问题，我在链接上 https://www.mathworks.com/matlab ... ty-is-close-to-zero 看到两种解决方法：一种是在MATLAB采用符号变量sym，一种是通过迭代的方法。这里不再赘述。下面说我的问题怎么解决，相应的MATLAB代码如下： Beta = u-；即是上面问题描述的截断的位置 % % method#1：经典方法 U1 = rand(1); X = norminv(U1+(1-U1).*normcdf(Beta)); % % method#2： 采用matlab的sym变量 U1 = rand(1); p = (U1+(1-U1).*(1-sym(normcdf(Beta,'upper')))); X = double(norminv(p)); % % method#3：文献链接https://arxiv.org/pdf/0907.4010.pdf，提供了一种采用平动指数分布函数作为抽样函数的Accept-Reject模拟方法，根据此编制了MATLAB函数 LeftTruncatedNormrnd X = LeftTruncatedNormrnd(Beta,Ns);， |
3楼2020-09-06 14:14:19
|
编制的MATLBA函数 LeftTruncatedNormrnd及其验证代码如下，经过验证感觉这种方法还不错。有路过大佬或虫友如果觉得有什么不对，或者更高效的方法，还请多指教。 主函数： function [X,accept_rate] = LeftTruncatedNormrnd(ul,Ns) alpha = (ul+sqrt(ul^2+4))/2; Ns_accept = 0; while Ns_accept==0 % 1. Generate translated exponential distribution Z = exprnd(1/alpha,Ns,1)+ul; % 2. Compute coeffcient of accept rate po = exp(-(Z-alpha).^2/2); % 3. Accept or reject U = rand(Ns,1); index = find(U<=po); X = Z(index,1); % 4. other Ns_accept = length(index); accept_rate = Ns_accept/Ns; end return 验证代码： xc_max = max(X); xc = ul  xc_max-ul)/1000:xc_max;dx = xc(2)-xc(1); n1 = hist(X,xc); pdf_sim = n1/Ns_accept/dx; pdf_norm = normpdf(xc,0,1)/normcdf(ul,'upper'); plot(xc,pdf_sim) hold on plot(xc,pdf_norm) |
4楼2020-09-06 14:18:34
5楼2020-09-06 14:20:39
|
验证代码更正，需先保存主函数LeftTruncatedNormrnd： ul = 2 Ns = 1e5; % method#1 % U1 = rand(Ns,1); % X = norminv(U1+(1-U1).*normcdf(ul)); % % % method#2 % U1 = rand(Ns,1); % p = (U1+(1-U1).*(1-sym(normcdf(ul,'upper')))); % X = double(norminv(p)); % % method#3 [X,accept_rate] = LeftTruncatedNormrnd(ul,Ns); close all xc_max = max(X); xc = ul xc_max-ul)/1000:xc_max;dx = xc(2)-xc(1); n1 = hist(X,xc); pdf_sim = n1/(Ns*accept_rate)/dx; pdf_norm = normpdf(xc,0,1)/normcdf(ul,'upper'); plot(xc,pdf_sim) hold on plot(xc,pdf_norm) legend('sim','target') |
6楼2020-09-06 14:48:54
