我正在使用 R 來嘗試求解
標準
正態分布函式
的公式。我遇到的問題是g的某些值非常小或在 R 中編碼為 0，因此給定的解決方案是 Inf。在下面的代碼行中給出了一個當g應該為正時編碼為 0 的示例。
pnorm(29.1,0,1) - pnorm(29,0,1) #this returns a value of 0
鑒于這個問題，我有兩個問題。首先，是否有一個數學技巧可以用來從上面的代碼中得到一個正數？我知道有一個pnorm引數log.p=TRUE可以回傳概率的對數。也許這會派上用場？
我的第二個問題是，在我有一個正但非常小的g的情況下，是否有一個數學技巧來確保f和g的比率的對數是有限的？
uj5u.com熱心網友回復：
正如答案和評論中指出的那樣，您實際上只需要計算log(g).
一種解決方案是使用對數空間加法/減法，它使用了一個很好的技巧來log(exp(x) ± exp(y))準確計算。
DPQ::logspace.sub(pnorm(29.1, log.p = TRUE), pnorm(29, log.p = TRUE))
## [1] -424.8435
@Alex 指出，如果我們通過反轉符號來查看（對稱）相反的尾部，以便我們從另一個非常小的數字中減去一個非常小的數字，則舍入誤差不是問題：
log(pnorm(-29,0,1) - pnorm(-29.1,0,1))
## -424.8435
比較來自@GregorThomas 的基于導數的估計：
dnorm(29, log = TRUE) log(0.1)
## [1] -423.7215
使用軟體包檢查這些解決方案Rmpfr：
library(Rmpfr)
x1 <- mpfr(29.1, 1000) ### 1000-bit precision
x2 <- mpfr(29.0, 1000)
log(pnorm(x1) - pnorm(x2))
## 1 'mpfr' number of precision 1000 bits
## [1] -424.843525917118620331088762368642820931865149932737039093516330881321230653835360169156688142607902967369067320263733953321675429119215305576014352760866238194504221709021981332048618426550583180796866417925580695487036011142823944123174363138216647098415472494704193913895461975221571930136768807682871
所以看起來@GregorThomas 的解決方案非常接近，但其他兩個解決方案更好。
uj5u.com熱心網友回復：
log(f / g) = log(f) - log(g)，所以你需要計算log(g)。
你有g = pnorm(x d) - pnorm(x).
密度函式 (dnorm()在 R 中) 是 CDF ( pnorm()) 的導數。這意味著dnorm(x) = (pnorm(x delta) - pnorm(x)) / delta隨著delta接近 0。所以對于小的delta我們可以近似dnorm(x) * delta = pnorm(x delta) - pnorm(x). 因此g(x) = dnorm(x) * d（大約）和log(g(x)) = log(dnorm(x)) log(d)。
在dnorm(29, log = TRUE) log(0.1)您的計算中用作log(g).
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/511536.html
標籤：r数学精确