[Stat] 數值 常態分佈累積機率函數 比較


最近開始整理水文統計網頁,打算重新用 ASP.NET 再寫一遍。 (關於 水文統計 – 線上分析 網頁)
最早的程式碼使用 Quick Basic 4.5 開發,是大四上 (1993) 修水文統計時開始寫的,當時以周文德[1]水文學原文為主。
念博一前半年 (2000) 的時候改寫成 ASP,採用 VBScript,那時工作要做動態網頁,所以當練功將大四的程式碼拿來改。

開始改寫時,想說精益求精,就開始把相關數值模型重新 research 。

參考維基百科[2] 的內容,原來數值模型還有好幾種,興起比較的意圖。

先用原先周文德水文學的方程式係數找到原作者跟年分[3],分別把方程式使用 Excel 2013 VBA 寫成函數比較,另外加入 Excel 內建的 Excel NORM.S.DIST 函數一起比較,如下表:

標準差 1 2 3 4 5 6
Wiki 表 0.8413447460685000 0.9772498680520000 0.9986501019685000 0.9999683287580000 0.9999997133485000 0.9999999990135000
Excel NORM.S.DIST 0.8413447460685430 0.9772498680518210 0.9986501019683700 0.9999683287581670 0.9999997133484280 0.9999999990134120
Abramowitz and Stegun, 1965
周文德水文學
0.8411238352705170 0.9774369986442540 0.9984208362961290 0.9999107486493130 0.9999941679910000 0.9999995054414810
Zelen and Severo, 1964 0.8413447404219760 0.9772499379836130 0.9986500327776380 0.9999683139654060 0.9999997128950000 0.9999999990098780
Hart, 1968 0.8413447460685430 0.9772498680518210 0.9986501019683700 0.9999683287581670 0.9999997133484280 0.9999999990134120
Marsaglia, 2004 0.8413447461005850 0.9772498680966180 0.9986500987983220 0.9999401630918050 0.9935899936665050 0.8832228594191580

由於不知道正確數值該多少,所以從維基百科[2]的表格做為比對依據,網頁上的數值位數有限,拿來跟 Excel 函數比對,假定 Excel 函數是正確的。

Marsaglia 使用泰勒級數展開,初步測試,17 項以上,小數點就沒有太大的變化,所以使用 20 項,2*20+1 = 41 ,x 的最高次方項為 x^41 。做到泰勒級數展開式時,不得不回來幹繳 VBA 7.1 LongLong 的識別字元是 ^ ,搞得我級數函數一值錯…

[VBA] LongLong 超長整數

Hart 則直接參考 West 2009 年的論文[4],裡面直接用 VB6 寫,剪來貼就可以了,West 說他把原始的 Fortran 函數轉成 VB 函數,果然我以前 Fortran 沒白學,以前的學者都用 Fortran 。

從比較表來看,以前水文統計使用周文德的累積機率函數誤差有點大,Zelen 的誤差都小於該函數,而 Hart 跟 Excel 結果完全一致,該不會 Excel 也用此函數來近似吧?

而 Marsaglia 使用泰勒級數展開,照理說泰勒級數展開應該要更準,不知道是寫錯還是其他問題,我懷疑也有可能是數值精度誤差問題,因為 x ^ 41 已經非常小,算完再去做除數,感覺誤差就被放大了。

所以我決定新版的水文統計在常態分佈累積機率函數就使用 Hart 了。(皮卡丘就是你了)

寫在這也希望給水利海洋相關從業人員參考,找個新函數比周文德的教科書更逼近,好吧,這是資訊從業人員的一種病。

1. https://zh.wikipedia.org/wiki/%E5%91%A8%E6%96%87%E5%BE%B7
2. https://en.wikipedia.org/wiki/Normal_distribution#Numerical_approximations_for_the_normal_CDF
3. https://books.google.com.tw/books?id=V1pHDwAAQBAJ&pg=PT523&lpg=PT523
4. https://s2.smu.edu/~aleskovs/emis/sqc2/accuratecumnorm.pdf

銀河的歷史,又翻過一頁~ 今年要重畫了,好期待~~~

廣告
Categories: 技術分享, 決策支援系統 | 標籤: | 發表留言

文章分頁導航

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com 標誌

您的留言將使用 WordPress.com 帳號。 登出 /  變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 /  變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 /  變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 /  變更 )

連結到 %s

在WordPress.com寫網誌.

%d 位部落客按了讚: