程式

Python中的三重高斯積分

  • November 29, 2021

我需要應用 scipy.stats.multivariate_normal.cdf(),它計算積分

$$ \int \frac{1}{\sqrt{(2\pi)^{\frac{3}{2}}\det(\Sigma)}} \exp\left(-\frac{1}{2}(x-\mu)^T \Sigma^{-1}(x-\mu)\right) dx $$,

在哪裡 $ \mu $ 是平均值並且 $ \Sigma^{-1} $ 是共變異數矩陣的逆矩陣 $ \Sigma $ .

根據我自己的理解,我試圖在另一個程序(Maple)中執行相同的計算。讓

\西格瑪 = $$ \begin{pmatrix} a_{11} & a_{12} & a_{13}\ a_{21} & a_{22} & a_{23}\ a_{31} & a_{32} & a_{33} \end{pmatrix} $$

為了這個計算的目的,讓 $ a_{11} = a_{22} = a_{33} = 1; a_{12} = a_{21} = 0.87055, a_{13} = a_{31} = 0.710804075, a_{23} = a_{32} = 0.8165 $ .

如果我們以指數形式擴展項,我們得到

$$ -\frac{1}{2}(x-\mu)^T \Sigma^{-1}(x-\mu) = -\frac{x^2}{2} - \frac{(y - 0.87055 x)^2}{2(1 - 0.87055^2)} - \frac{(z - 0.8165 y)^2}{2(1 - 0.8165^2)} $$

所以我試圖在 Maple 中評估的積分是:

$$ \int\limits_{-\infty}^{1.37824} \int\limits_{-\infty}^{-21.58961} \int\limits_{-\infty}^{18.48617} \frac{\exp\left( -\frac{x^2}{2} - \frac{(y - 0.87055 x)^2}{2(1 - 0.87055^2)} - \frac{(z - 0.8165 y)^2}{2(1 - 0.8165^2)} \right)}{(2\pi)^{\frac{3}{2}}\sqrt{(1-0.87055^2)(1-0.8165^2)}} dz dy dx $$.

結果是 $ 4.164024864\times10^{-242} $ .

應用 scipy.stats.multivariate_normal.cdf() 與 $ [x, y, z] = [1.37824, -21.58961, 18.48617] $ , $ mean=None $ 和 $ cov=np.array([[1, 0.87055, 0.710804075], [0.87055, 1, 0.8165], [0.710804075, 0.8165, 1]]) $ .

Python中的結果是 $ 1.124512788731174 \times 10^{-103} $ .

兩個結果之間的差異是顯著的。顯然我做錯了什麼。如果有人能指出我的錯誤在哪裡,我將不勝感激!先感謝您!

這兩個值實際上都為零,並且兩個軟體之間的差異非常非常小(與您的輸入相比幾乎是無窮小)。電腦不是完美的計算器,並且會出現一些舍入誤差。嘗試注意它們可能發生並檢查舍入誤差是否對您的應用程序很重要。

引用自:https://quant.stackexchange.com/questions/68928