請(qǐng)教如何提高M(jìn)ATLAB函數(shù)norminv(p,0,1)的使用精度
請(qǐng)教各位大佬一個(gè)MATLAB問(wèn)題：如何提高函數(shù)norminv(p,0,1)的使用精度？
norminv為一維正態(tài)分布的累積分布函數(shù)的逆函數(shù)，假設(shè) p = 1-10^(i) , i = 1,...,.N， 因?yàn)閙atlab的運(yùn)算精度為雙精度16位有效數(shù)字 ，因此當(dāng)N>16的時(shí)候，norminv(p,0,1)返回的都是正無(wú)窮大�，F(xiàn)在編程涉及到的p比較小，所以想請(qǐng)教一下如何使得p<1-10^16時(shí)，也能使得該函數(shù)得出較為準(zhǔn)確的結(jié)果？
返回小木蟲(chóng)查看更多
今日熱帖
京公網(wǎng)安備 11010802022153號(hào)
做數(shù)值計(jì)算的第一步就應(yīng)該是把問(wèn)題歸一化，從而在計(jì)算質(zhì)量最高的范圍求解
該問(wèn)題已經(jīng)解決。實(shí)際上，我本來(lái)的目的是為生成服從左截?cái)嗟恼龖B(tài)分布變量，X~N(mu,sigma,u-), 式中mu,sigma為標(biāo)準(zhǔn)正態(tài)變量的均值與方差，u-為左截?cái)嗟奈恢�。通常的的方法是生成均勻分布變�?p~U(p1,1)，其中p1 = normcdf(u-,mu,sigma)，然后再根據(jù)正態(tài)函數(shù)的累積分布函數(shù)將p1轉(zhuǎn)化目標(biāo)的截?cái)嗾龖B(tài)分布變量X， 即是X = norminv(p1,mu,sigma)。其中normcdf，norminv都是MATLAB內(nèi)置函數(shù)，分別是正態(tài)分布函數(shù)的累積分布函數(shù)與其逆函數(shù)。這種做法在u-與mu比較接近的時(shí)候，效率還是可以的。但是當(dāng)u->>mu，即是p1非常接近于1，如同我在帖子中所描述的問(wèn)題，由于數(shù)值計(jì)算的精度所限制這種做法是行不通的。
關(guān)于當(dāng)p1非常接近于1時(shí)，對(duì)應(yīng)的逆函數(shù)值的求解問(wèn)題，我在鏈接上 https://www.mathworks.com/matlab ... ty-is-close-to-zero 看到兩種解決方法：一種是在MATLAB采用符號(hào)變量sym，一種是通過(guò)迭代的方法。這里不再贅述。下面說(shuō)我的問(wèn)題怎么解決，相應(yīng)的MATLAB代碼如下：
Beta = u-；即是上面問(wèn)題描述的截?cái)嗟奈恢?br /> % % method#1：經(jīng)典方法
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：文獻(xiàn)鏈接https://arxiv.org/pdf/0907.4010.pdf，提供了一種采用平動(dòng)指數(shù)分布函數(shù)作為抽樣函數(shù)的Accept-Reject模擬方法，根據(jù)此編制了MATLAB函數(shù) LeftTruncatedNormrnd
X = LeftTruncatedNormrnd(Beta,Ns);，
編制的MATLBA函數(shù) LeftTruncatedNormrnd及其驗(yàn)證代碼如下，經(jīng)過(guò)驗(yàn)證感覺(jué)這種方法還不錯(cuò)。有路過(guò)大佬或蟲(chóng)友如果覺(jué)得有什么不對(duì)，或者更高效的方法，還請(qǐng)多指教。
主函數(shù)：
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
驗(yàn)證代碼：
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)
代碼中哭臉：( 為 ‘:’ + '('，請(qǐng)自行替換，比較少發(fā)帖，還請(qǐng)見(jiàn)諒
驗(yàn)證代碼更正，需先保存主函數(shù)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')