Saya telah melihat simulasi Monte Carlo baru-baru ini, dan telah menggunakannya untuk memperkirakan konstanta seperti (lingkaran di dalam persegi panjang, area proporsional).$\pi$

Namun, saya tidak dapat memikirkan metode yang sesuai untuk memperkirakan nilai [angka Euler] menggunakan integrasi Monte Carlo.$e$

Apakah Anda memiliki petunjuk tentang bagaimana hal ini dapat dilakukan?

Cara sederhana dan elegan untuk memperkirakan oleh Monte Carlo dijelaskan dalam makalah ini . Makalah ini sebenarnya tentang pengajaran . Karenanya, pendekatan tersebut tampaknya sangat sesuai untuk tujuan Anda. Idenya didasarkan pada latihan dari buku teks Rusia populer tentang teori probabilitas oleh Gnedenko. Lihat contoh.22 di hlm.183$e$$e$

Itu terjadi sehingga , di mana adalah variabel acak yang didefinisikan sebagai berikut. Ini adalah jumlah minimum sehingga dan adalah angka acak dari distribusi seragam pada . Cantik bukan ?!ξ n ∑ n i = 1 r i > 1 r i [ 0 , 1 $E[\xi]=e$$\xi$$n$$\sum_{i=1}^n r_i>1$$r_i$$[0,1]$

Karena ini adalah latihan, saya tidak yakin apakah itu keren bagi saya untuk mengirim solusi (bukti) di sini :) Jika Anda ingin membuktikannya sendiri, inilah tipnya: bab ini disebut "Momen", yang seharusnya menunjukkan Anda ke arah yang benar.

Jika Anda ingin menerapkannya sendiri, maka jangan membaca lebih lanjut!

Ini adalah algoritma sederhana untuk simulasi Monte Carlo. Gambarlah seragam acak, lalu yang lain dan seterusnya sampai jumlahnya melebihi 1. Jumlah tebusan yang ditarik adalah percobaan pertama Anda. Katakanlah Anda mendapat:

 0.0180
 0.4596
 0.7920

Kemudian sidang pertama Anda diberikan 3. Jauhkan melakukan percobaan ini, dan Anda akan melihat bahwa rata-rata Anda mendapatkan .$e$

Kode MATLAB, hasil simulasi dan histogram mengikuti.

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)

Hasil dan histogram:

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

UPDATE: Saya memperbarui kode saya untuk menyingkirkan berbagai hasil uji coba sehingga tidak memerlukan RAM. Saya juga mencetak estimasi PMF.

Pembaruan 2: Inilah solusi Excel saya. Masukkan tombol di Excel dan tautkan ke makro VBA berikut:

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

Masukkan jumlah uji coba, seperti 1000, di sel D1, dan klik tombol. Di sini bagaimana tampilan layar setelah menjalankan pertama:

UPDATE 3: Silverfish menginspirasi saya ke cara lain, tidak seanggun yang pertama tapi masih keren. Ini menghitung volume n-simpleks menggunakan urutan Sobol .

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

Secara kebetulan dia menulis buku pertama tentang metode Monte Carlo yang saya baca di sekolah menengah. Ini pengantar metode terbaik menurut saya.

PEMBARUAN 4:

Silverfish dalam komentar menyarankan implementasi rumus Excel sederhana. Ini adalah jenis hasil yang Anda dapatkan dengan pendekatannya setelah sekitar 1 juta angka acak dan uji coba 185 ribu:

Jelas, ini jauh lebih lambat daripada implementasi Excel VBA. Terutama, jika Anda mengubah kode VBA saya untuk tidak memperbarui nilai sel di dalam loop, dan hanya lakukan setelah semua statistik dikumpulkan.

PEMBARUAN 5

Solusi Xi'an # 3 terkait erat (atau bahkan sama dalam arti sesuai komentar jwg di utas). Sulit untuk mengatakan siapa yang datang dengan gagasan pertama Forsythe atau Gnedenko. Edisi 1950 asli Gnedenko dalam bahasa Rusia tidak memiliki bagian Masalah di Bab. Jadi, saya tidak dapat menemukan masalah ini pada pandangan pertama di mana ia berada di edisi selanjutnya. Mungkin ditambahkan kemudian atau dikubur dalam teks.

Seperti yang saya komentari dalam jawaban Xi'an, pendekatan Forsythe terkait dengan bidang menarik lainnya: distribusi jarak antara puncak (ekstrem) dalam urutan acak (IID). Jarak rata-rata adalah 3. Urutan ke bawah dalam pendekatan Forsythe berakhir dengan dasar, jadi jika Anda melanjutkan pengambilan sampel Anda akan mendapatkan bagian bawah lainnya di beberapa titik, lalu di titik lainnya. Anda dapat melacak jarak antara mereka dan membangun distribusi.

Saya menyarankan untuk tidak memilih jawaban Aksakal. Ini tidak bias dan hanya bergantung pada metode menghasilkan penyimpangan seragam unit.

Jawaban saya dapat dibuat tepat secara sewenang-wenang, tetapi masih bias dari nilai sebenarnya dari .$e$

Jawaban Xi'an benar, tetapi saya pikir ketergantungannya pada fungsi atau cara menghasilkan penyimpangan Poisson acak agak melingkar ketika tujuannya adalah untuk memperkirakan e .$\log$$e$

Memperkirakan dengan Bootstrapping$e$

Sebagai gantinya, pertimbangkan prosedur bootstrap. Satu memiliki sejumlah besar objek yang ditarik dengan penggantian ke ukuran sampel n . Pada setiap undian, probabilitas untuk tidak menggambar objek i tertentu adalah 1 - n - 1 , dan ada n gambar seperti itu. Probabilitas bahwa objek tertentu dihilangkan dari semua undian adalah p = ( 1 - $n$$n$$i$$1-n^{-1}$$n$$p=(1-\frac{1}{n})^n.$

jadi kita juga dapat menulis

Yaitu, estimasi ditemukan dengan memperkirakan probabilitas bahwa pengamatan khusus dihilangkan dari bootstrap mereplikasi melintasi banyak ulangan seperti itu - yaitu fraksi kemunculan objek dalam bootstraps.m B $p$$m$$B_j$$i$

Ada dua sumber kesalahan dalam perkiraan ini. Hingga akan selalu berarti bahwa hasilnya adalah perkiraan, yaitu estimasi bias. Selain itu, akan berfluktuasi di sekitar nilai sebenarnya karena ini adalah simulasi.$n$$\hat{p}$

Saya menemukan pendekatan ini agak menarik karena sarjana atau orang lain dengan cukup sedikit untuk melakukan bisa perkiraan menggunakan setumpuk kartu, tumpukan batu-batu kecil, atau item lainnya di tangan, di vena sama seperti seseorang bisa memperkirakan menggunakan kompas, tepi lurus dan beberapa butir pasir. Saya pikir itu rapi ketika matematika dapat dipisahkan dari kenyamanan modern seperti komputer.$e$$\pi$

Hasil

Saya melakukan beberapa simulasi untuk berbagai jumlah replikasi bootstrap. Kesalahan standar diperkirakan menggunakan interval normal.

Perhatikan bahwa pilihan jumlah objek yang sedang bootstrap menetapkan batas atas mutlak pada keakuratan hasil karena prosedur Monte Carlo memperkirakan dan hanya bergantung pada . Pengaturan menjadi terlalu besar hanya akan membebani komputer Anda, baik karena Anda hanya memerlukan perkiraan "kasar" ke atau karena bias akan dibanjiri oleh varians karena Monte Carlo. Hasil ini untuk dan akurat untuk desimal ketiga.p p n n e n = 10 3 p - 1 ≈ $n$$p$$p$$n$$n$$e$$n=10^3$$p^{-1}\approx e$

Plot ini menunjukkan bahwa pilihan memiliki konsekuensi langsung dan mendalam bagi stabilitas di . Garis putus-putus biru menunjukkan dan garis merah menunjukkan . Seperti yang diharapkan, meningkatkan ukuran sampel menghasilkan perkiraan yang lebih akurat . p p e $m$$\hat{p}$$p$$e$$\hat{p}$

Saya menulis naskah R yang panjang dan memalukan untuk ini. Saran untuk perbaikan dapat diajukan di belakang tagihan $ 20.

library(boot)
library(plotrix)
n <- 1e3

## if p_hat is estimated with 0 variance (in the limit of infinite bootstraps), then the best estimate we can come up with is biased by exactly this much:
approx <- 1/((1-1/n)^n)

dat <- c("A", rep("B", n-1))
indicator <- function(x, ndx)   xor("A"%in%x[ndx], TRUE) ## Because we want to count when "A" is *not* in the bootstrap sample

p_hat <- function(dat, m=1e3){
    foo <- boot(data=dat, statistic=indicator, R=m) 
    1/mean(foo$t)
} 

reps <- replicate(100, p_hat(dat))

boxplot(reps)
abline(h=exp(1),col="red")

p_mean <- NULL
p_var <- NULL
for(i in 1:10){
    reps <- replicate(2^i, p_hat(dat))
    p_mean[i] <- mean(reps)
    p_var[i] <- sd(reps)
}
plotCI(2^(1:10), p_mean, uiw=qnorm(0.975)*p_var/sqrt(2^(1:10)),xlab="m", log="x", ylab=expression(hat(p)), main=expression(paste("Monte Carlo Estimates of ", tilde(e))))
abline(h=approx, col='red')

Solusi 1:

Untuk distribusi Poisson , Oleh karena itu, jika , yang berarti Anda dapat memperkirakan oleh simulasi Poisson. Dan simulasi Poisson dapat diturunkan dari generator distribusi eksponensial (jika tidak dengan cara yang paling efisien).P ( X = k ) = λ $\mathcal{P}(\lambda)$ X ∼ P ( 1 ) P ( X = 0 ) = P ( X = 1 ) = e - 1 e -

Catatan 1: Sebagaimana dibahas dalam komentar, ini adalah argumen yang agak berbelit-belit karena mensimulasikan dari distribusi Poisson atau ekuivalen distribusi Eksponensial mungkin sulit dibayangkan tanpa melibatkan log atau fungsi exp ... Tapi kemudian W. Huber datang ke selamatkan jawaban ini dengan solusi paling elegan berdasarkan seragam yang dipesan. Namun yang merupakan perkiraan , karena distribusi spasi seragam adalah Beta , menyiratkan bahwa yang konvergen ke sebagai$U_{(i:n)}-U_{(i-1:n)}$$\mathfrak{B}(1,n)$

$$\mathbb{P}(n\{U_{(i:n)}-U_{(i-1:n)}\}\ge 1)=\left(1-\frac{1}{n}\right)^n$$ $e^{-1}$$n$tumbuh hingga tak terbatas. Sebagai samping lain yang menjawab komentar, generator eksponensial tahun 1951 von Neumann hanya menggunakan generasi yang seragam.

Solusi 2:

Cara lain untuk mencapai representasi konstanta sebagai integral adalah dengan mengingat bahwa, ketika lalu yang juga merupakan distribusi . Oleh karena itu, Pendekatan kedua untuk mendekati oleh Monte Carlo dengan demikian mensimulasikan pasangan normal dan memantau frekuensi kali . Dalam arti itu adalah kebalikan dari perkiraan Monte Carlo dari terkait dengan frekuensi kali ...$e$

Solusi 3:

Rekan saya dari Universitas Warwick, M. Pollock, menunjukkan perkiraan lain dari Monte Carlo yang disebut metode Forsythe : idenya adalah menjalankan serangkaian generasi yang seragam sampai . Ekspektasi aturan pemberhentian yang sesuai, , yang merupakan jumlah waktu urutan seragam turun maka sedangkan probabilitas bahwa ganjil adalah ! ( Metode Forsythe sebenarnya bertujuan mensimulasikan dari setiap kepadatan bentuk , maka lebih umum daripada mendekati dan .)$u_1,u_2,...$$u_{n+1}>u_{n}$$N$$e$$N$$e^{-1}$$\exp G(x)$$e$$e^{-1}$

Ini sangat paralel dengan pendekatan Gnedenko yang digunakan dalam jawaban Aksakal , jadi saya ingin tahu apakah yang satu dapat berasal dari yang lain. Paling tidak, keduanya memiliki distribusi yang sama dengan massa probabilitasuntuk nilai .$1/n!$$n$

Implementasi R cepat dari metode Forsythe adalah untuk tidak mengikuti dengan tepat urutan seragam yang mendukung blok yang lebih besar, yang memungkinkan untuk pemrosesan paralel:

use=runif(n)
band=max(diff((1:(n-1))[diff(use)>0]))+1
bends=apply(apply((apply(matrix(use[1:((n%/%band)*band)],nrow=band),
2,diff)<0),2,cumprod),2,sum)

Bukan solusi ... hanya komentar cepat yang terlalu panjang untuk kotak komentar.

Aksakal

Aksakal memposting solusi di mana kami menghitung jumlah yang diharapkan dari gambar Seragam standar yang harus diambil, sehingga jumlahnya akan melebihi 1. Dalam Mathematica , formulasi pertama saya adalah:

mrM := NestWhileList[(Random[] + #) &, Random[], #<1 &]

Mean[Table[Length[mrM], {10^6}]] 

EDIT: Baru saja bermain cepat dengan ini, dan kode berikut (metode yang sama - juga di Mma - kode yang berbeda) sekitar 10 kali lebih cepat:

Mean[Table[Module[{u=Random[], t=1},  While[u<1, u=Random[]+u; t++]; t] , {10^6}]]

Xian / Whuber

Whuber telah menyarankan kode keren cepat untuk mensimulasikan solusi Xian 1:

Versi R: n <- 1e5; 1/mean(n*diff(sort(runif(n+1))) > 1)

Versi mma: n=10^6; 1. / Mean[UnitStep[Differences[Sort[RandomReal[{0, n}, n + 1]]] - 1]]

yang dia catat adalah 20 kali lebih cepat kode pertama (atau sekitar dua kali lebih cepat dari kode baru di atas).

Hanya untuk bersenang-senang, saya pikir akan menarik untuk melihat apakah kedua pendekatan sama efisiennya (dalam arti statistik). Untuk melakukannya, saya menghasilkan 2000 perkiraan e menggunakan:

• Metode Aksakal: dataA
• Metode 1 Xian menggunakan kode whuber: dataB

... keduanya dalam Mathematica . Diagram berikut kontras estimasi kerapatan kernel nonparametrik dari data yang dihasilkan dan dataB set data.

Jadi, sementara kode whuber (kurva merah) sekitar dua kali lebih cepat, metode ini tampaknya tidak dapat diandalkan.

Metode yang membutuhkan jumlah sampel yang tidak saleh

$f(x) = e^x$$\frac{ \bar x \sqrt{12}}{ \sqrt{n}} \dot \sim N(0,1)$$e$$\infty$$\infty$

$N(0,1)$$e^{x}$

$N(0,1)$$\hat \phi(x)$$\phi(\sqrt(2) ) = (2 \pi)^{-1/2} e^{-1}$$e = \hat \phi(\sqrt{2}) \sqrt{2 \pi}$

Jika Anda ingin benar-benar gila, Anda bahkan dapat memperkirakan dan menggunakan metode yang Anda bahas sebelumnya.$\sqrt{2}$$\sqrt{2\pi}$

Metode yang membutuhkan sampel sangat sedikit, tetapi menyebabkan jumlah kesalahan numerik yang tidak baik

Jawaban yang benar-benar konyol, tetapi sangat efisien, berdasarkan komentar yang saya buat:

Biarkan . Tentukan. Definisikan .Y n = | ( ˉ x ) n | E = ( 1 - Y n ) - 1 / Y $X \sim \text{uniform}(-1, 1)$$Y_n = | (\bar x)^n|$$\hat e = (1 - Y_n)^{-1/Y_n}$

Ini akan menyatu sangat cepat, tetapi juga mengalami kesalahan numerik ekstrem.

1 / Y n n → ∞ Y n Y n = $Y_n$$1/Y_n$$n \rightarrow \infty$$Y_n$$Y_n = 0$$e$

Berikut cara lain yang bisa dilakukan, meskipun cukup lambat. Saya tidak mengklaim efisiensi, tetapi menawarkan alternatif ini dalam semangat kelengkapan.

$n$$U_1, \cdots , U_n \sim \text{IID U}(0,1)$$e$$^\dagger$

$e$$u_{(1)} \geqslant \cdots \geqslant u_{(n)}$

$m \equiv \min \{ k | S(k) \geqslant 1 \}$$1/e$$e$

$1/e$$e$$e$

Implementasi dalam R: Metode ini dapat diterapkan dalam Rmenggunakan runifuntuk menghasilkan nilai yang seragam. Kode tersebut adalah sebagai berikut:

EST_EULER <- function(n) { U <- sort(runif(n), decreasing = TRUE);
                           S <- cumsum(1/U)/n;
                           m <- min(which(S >= 1));
                           2/(U[m-1]+U[m]); }

$e$

set.seed(1234);

EST_EULER(10^3);
[1] 2.715426

EST_EULER(10^4);
[1] 2.678373

EST_EULER(10^5);
[1] 2.722868

EST_EULER(10^6); 
[1] 2.722207

EST_EULER(10^7);
[1] 2.718775

EST_EULER(10^8);
[1] 2.718434

> exp(1)
[1] 2.718282

$^\dagger$$e$