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

Euler-Mascheroni 常數被簡單地定義為調和級數與自然對數之間的極限差。 γ=limn→∞(n∑k=11k−lnn)

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

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

我認為使用Gumbel Distribution可能是可行的 μ=0,β=1 （因為 Gumbel 分佈的平均值是 E(X)=μ+βγ ）。但是，我不確定如何實現這一點。

**編輯：**感謝 Xi’an 和 S. Catterall 指出模擬 Gumbel Dist 的說明在 Wikipedia 文章本身上。 Q(p)=μ−βln(−ln(p)),

你畫的地方 p 從 (0,1) .

簡單的 Mathematica 代碼 106 畫，

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

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

γ=−∞∫0e−xlogx dx\[6pt]=∞∫0(−logx)⋅pX(x) dx\[12pt]=E(−logX).\[12pt]

因此，我們可以採取 X1,…,Xn∼IID Exp(1) 並使用估算器：

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

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

#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