Simulation

估計 Euler-Mascheroni 常數 (γγgamma) 通過蒙特卡洛模擬

  • June 21, 2021

Euler-Mascheroni 常數被簡單地定義為調和級數與自然對數之間的極限差。 $$ \gamma =\lim_{n\to \infty}\left(\sum _{k=1}^{n}{\frac {1}{k}}-\ln n\right) $$

在看到以下帖子Approximate後,我對使用 Monte Carlo 模擬的數學常數的估計很感興趣 使用蒙特卡洛模擬並閱讀無數實驗來近似 $ \pi $ (布馮的針/麵條等)。我想展示一些 MC 模擬計算數學常數的好例子,以幫助激發這個想法。

我的問題是, 你如何設計一個非平凡的 MC 模擬來估計 $ \gamma $ ? (我說很重要,因為我可以整合 $ \ln n $ 對於足夠大的 $ n $ 使用數字 MC 積分,但這似乎不是一個非常有趣的展示。)

我認為使用Gumbel Distribution可能是可行的 $ \mu=0, \beta=1 $ (因為 Gumbel 分佈的平均值是 $ E(X)=\mu+\beta\gamma $ )。但是,我不確定如何實現這一點。

**編輯:**感謝 Xi’an 和 S. Catterall 指出模擬 Gumbel Dist 的說明在 Wikipedia 文章本身上。 $$ {\displaystyle Q(p)=\mu -\beta \ln(-\ln(p)),} $$你畫的地方 $ p $ 從 $ (0,1) $ .

簡單的 Mathematica 代碼 $ 10^6 $ 畫,

gumbelrand = RandomReal[{0, 1}, {10^6, 1}];
Mean[-Log[-Log[gumbelrand]]]

如果您願意為您的方法使用指數和對數函數,您可以 使用指數隨機變量的重要性抽樣來估計Euler-Mascheroni 常數。讓 $ X \sim \text{Exp}(1) $ 我們可以將 Euler-Mascheroni 常數寫為:

$$ \begin{align} \gamma &= - \int \limits_0^\infty e^{-x} \log x \ dx \[6pt] &= \int \limits_0^\infty (-\log x) \cdot p_X(x) \ dx \[12pt] &= \mathbb{E}(- \log X). \[12pt] \end{align} $$

因此,我們可以採取 $ X_1,…,X_n \sim \text{IID Exp}(1) $ 並使用估算器:

$$ \hat{\gamma}n \equiv - \frac{1}{n} \sum{i=1}^n \log x_i. $$

大數定律確保這是一個強一致的估計量 $ \gamma $ , 所以它會隨機收斂到這個值 $ n \rightarrow \infty $ . 我們可以R很簡單地實現這一點,如下所示。使用 $ n=10^6 $ 模擬我們得到非常接近真實值。

#Create function to estimate Euler–Mascheroni constant
simulate.EM <- function(n) { -mean(log(rexp(n, rate = 1))) }

#Perform simulations
set.seed(1)
simulate.EM(10^6)
[1] 0.5772535

#True value is 0.5772157
#Our simulation is correct to four DP

引用自:https://stats.stackexchange.com/questions/531538

comments powered by Disqus