请教如何提高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时，也能使得该函数得出较为准确的结果？
返回小木虫查看更多
今日热帖
京公网安备 11010802022153号
做数值计算的第一步就应该是把问题归一化，从而在计算质量最高的范围求解
该问题已经解决。实际上，我本来的目的是为生成服从左截断的正态分布变量，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);，
编制的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 = ulxc_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)
代码中哭脸：( 为 ‘:’ + '('，请自行替换，比较少发帖，还请见谅
验证代码更正，需先保存主函数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 = ulxc_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')