| 24小時(shí)熱門(mén)版塊排行榜 |
| 5 | 1/1 | 返回列表 |
| 查看: 2599 | 回復(fù): 5 | |||
| 【懸賞金幣】回答本帖問(wèn)題,作者dsmj1995將贈(zèng)送您 10 個(gè)金幣 | |||
| 當(dāng)前只顯示滿(mǎn)足指定條件的回帖,點(diǎn)擊這里查看本話(huà)題的所有回帖 | |||
[求助]
請(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ú)窮大。現(xiàn)在編程涉及到的p比較小,所以想請(qǐng)教一下如何使得p<1-10^16時(shí),也能使得該函數(shù)得出較為準(zhǔn)確的結(jié)果? |
|
編制的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 = 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) |
木蟲(chóng) (著名寫(xiě)手)
|
做數(shù)值計(jì)算的第一步就應(yīng)該是把問(wèn)題歸一化,從而在計(jì)算質(zhì)量最高的范圍求解 發(fā)自小木蟲(chóng)IOS客戶(hù)端 |

|
該問(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);, |
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 材料專(zhuān)碩306英一數(shù)二 +7 | z1z2z3879 2026-03-16 | 9/450 |
|
|---|---|---|---|---|
|
[考研] 344求調(diào)劑 +4 | knight344 2026-03-16 | 4/200 |
|
|
[考研] 本人考085602 化學(xué)工程 專(zhuān)碩 +16 | 不知道叫什么! 2026-03-15 | 18/900 |
|
|
[考研] 332求調(diào)劑 +6 | Zz版 2026-03-13 | 6/300 |
|
|
[考研] 277調(diào)劑 +4 | 自由煎餅果子 2026-03-16 | 5/250 |
|
|
[考研] 有沒(méi)有道鐵/土木的想調(diào)劑南林,給自己招師弟中~ +3 | TqlXswl 2026-03-16 | 7/350 |
|
|
[考研] 274求調(diào)劑 +5 | 時(shí)間點(diǎn) 2026-03-13 | 5/250 |
|
|
[考研]
|
zhouzhen654 2026-03-16 | 3/150 |
|
|
[文學(xué)芳草園] 伙伴們,祝我生日快樂(lè)吧 +17 | myrtle 2026-03-10 | 26/1300 |
|
|
[基金申請(qǐng)]
今年的國(guó)基金是打分制嗎?
50+3
|
zhanghaozhu 2026-03-14 | 3/150 |
|
|
[考研] 070300化學(xué)學(xué)碩求調(diào)劑 +6 | 太想進(jìn)步了0608 2026-03-16 | 6/300 |
|
|
[考研] 277材料科學(xué)與工程080500求調(diào)劑 +3 | 自由煎餅果子 2026-03-16 | 3/150 |
|
|
[教師之家] 焦慮 +7 | 水冰月月野兔 2026-03-13 | 9/450 |
|
|
[考博] 東華理工大學(xué)化材專(zhuān)業(yè)26屆碩士博士申請(qǐng) +6 | zlingli 2026-03-13 | 6/300 |
|
|
[考研] 22408總分284求調(diào)劑 +3 | InAspic 2026-03-13 | 3/150 |
|
|
[考研] 288求調(diào)劑 +4 | 奇點(diǎn)0314 2026-03-14 | 4/200 |
|
|
[考研] 0856材料與化工301求調(diào)劑 +5 | 奕束光 2026-03-13 | 5/250 |
|
|
[考研] 工科278分求調(diào)劑 +5 | 周慢熱啊 2026-03-12 | 7/350 |
|
|
[考研] 一志愿山大07化學(xué) 332分 四六級(jí)已過(guò) 本科山東雙非 求調(diào)劑! +3 | 不想理你 2026-03-12 | 3/150 |
|
|
[考研] 0856化學(xué)工程280分求調(diào)劑 +4 | shenzxsn 2026-03-11 | 4/200 |
|