近似𝑒ee使用蒙特卡羅模擬

February 4, 2016

我最近一直在研究蒙特卡洛模擬，並一直在使用它來近似常數，例如（矩形內的圓圈，比例區域）。

但是，我想不出相應的方法來近似[歐拉數] 使用蒙特卡洛積分。

您對如何做到這一點有任何指示嗎？

簡單而優雅的估算方法 $ e $ Monte Carlo 在本文中進行了描述。這篇論文實際上是關於教學的 $ e $ . 因此，該方法似乎非常適合您的目標。這個想法是基於 Gnedenko 的一本流行的烏克蘭概率論教科書的練習。見第 183 頁的 ex.22

它的發生是這樣的 $ E[\xi]=e $ ，在哪裡 $ \xi $ 是一個隨機變量，定義如下。這是最小數量 $ n $ 這樣 $ \sum_{i=1}^n r_i>1 $ 和 $ r_i $ 是來自均勻分佈的隨機數 $ [0,1] $ . 漂亮，不是嗎？！

由於這是一個練習，我不確定我在這裡發布解決方案（證明）是否很酷：）如果您想自己證明，這裡有一個提示：本章稱為“時刻”，應該指出你在正確的方向。

如果您想自己實現它，請不要進一步閱讀！

這是一個簡單的蒙特卡羅模擬算法。抽一個統一的隨機數，然後再抽一個，依此類推，直到總和超過 1。抽出的隨機數是您的第一次試驗。假設你有：
0.0180
0.4596
0.7920
然後你的第一個試驗呈現 3. 繼續做這些試驗，你會注意到平均你得到 $ e $ .

MATLAB 代碼、仿真結果和直方圖如下。
N = 10000000;
n = N;
s = 0;
i = 0;
maxl = 0;
f = 0;
while n > 0
   s = s + rand;
   i = i + 1;
   if s > 1
       if i > maxl
           f(i) = 1;
           maxl = i;
       else
           f(i) = f(i) + 1;
       end
       i = 0;
       s = 0;
       n = n - 1;
   end
end

disp ((1:maxl)*f'/sum(f))
bar(f/sum(f))
grid on

f/sum(f)
結果和直方圖：
2.7183


ans =

 Columns 1 through 8

        0    0.5000    0.3332    0.1250    0.0334    0.0070    0.0012    0.0002

 Columns 9 through 11

   0.0000    0.0000    0.0000
更新：我更新了我的代碼以擺脫試驗結果數組，這樣它就不會佔用 RAM。我還打印了 PMF 估計值。

更新 2：這是我的 Excel 解決方案。在 Excel 中放置一個按鈕並將其鏈接到以下 VBA 宏：
Private Sub CommandButton1_Click()
n = Cells(1, 4).Value
Range("A:B").Value = ""
n = n
s = 0
i = 0
maxl = 0
Cells(1, 2).Value = "Frequency"
Cells(1, 1).Value = "n"
Cells(1, 3).Value = "# of trials"
Cells(2, 3).Value = "simulated e"
While n > 0
   s = s + Rnd()
   i = i + 1
   If s > 1 Then
       If i > maxl Then
           Cells(i, 1).Value = i
           Cells(i, 2).Value = 1
           maxl = i
       Else
           Cells(i, 1).Value = i
           Cells(i, 2).Value = Cells(i, 2).Value + 1
       End If
       i = 0
       s = 0
       n = n - 1
   End If
Wend


s = 0
For i = 2 To maxl
   s = s + Cells(i, 1) * Cells(i, 2)
Next


Cells(2, 4).Value = s / Cells(1, 4).Value

Rem bar (f / Sum(f))
Rem grid on

Rem f/sum(f)

End Sub
在單元格 D1 中輸入試驗次數，例如 1000，然後單擊按鈕。這是第一次運行後屏幕的樣子：

更新 3： Silverfish 以另一種方式啟發了我，雖然不像第一個那樣優雅，但仍然很酷。它使用Sobol序列計算了 n-單純形的體積。
s = 2;
for i=2:10
   p=sobolset(i);
   N = 10000;
   X=net(p,N)';
   s = s + (sum(sum(X)<1)/N);
end
disp(s)

2.712800000000001
巧合的是，他寫了我在高中讀過的第一本關於蒙特卡洛方法的書。在我看來，這是對該方法的最佳介紹。

更新 4：

評論中的 Silverfish 建議了一個簡單的 Excel 公式實現。這是你在大約 100 萬個隨機數和 185K 次試驗後使用他的方法得到的結果：

顯然，這比 Excel VBA 實現要慢得多。特別是，如果您修改我的 VBA 代碼以不更新循環內的單元格值，並且僅在收集所有統計信息後才執行此操作。

更新 5

西安的解決方案#3 密切相關（甚至在某種意義上與 jwg 在線程中的評論相同）。很難說是誰首先提出了這個想法，Forsythe 或 Gnedenko。Gnedenko 的 1950 年俄文原版在章節中沒有問題部分。所以，乍一看，我在以後的版本中找不到這個問題。也許它是後來添加的或埋在文本中。

正如我在西安的回答中評論的那樣，Forsythe 的方法與另一個有趣的領域有關：隨機 (IID) 序列中峰值（極值）之間的距離分佈。平均距離恰好是 3。Forsythe 方法中的向下序列以底部結束，因此如果您繼續採樣，您將在某個點得到另一個底部，然後是另一個底部，等等。您可以跟踪它們之間的距離並建立分佈。