估計 Euler-Mascheroni 常數 (γγgamma) 通過蒙特卡洛模擬
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