Ostatnio przyglądałem się symulacji Monte Carlo i używałem jej do przybliżania stałych, takich jak (okrąg wewnątrz prostokąta, obszar proporcjonalny).$\pi$

Nie jestem jednak w stanie wymyślić odpowiedniej metody aproksymacji wartości $e$ [liczby Eulera] przy użyciu integracji Monte Carlo.

Czy masz jakieś wskazówki, jak to zrobić?

W tym artykule opisano prosty i elegancki sposób oszacowania przez Monte Carlo . Artykuł dotyczy nauczania . Dlatego podejście wydaje się idealnie pasować do twojego celu. Pomysł opiera się na ćwiczeniu z popularnego rosyjskiego podręcznika teorii prawdopodobieństwa autorstwa Gnedenko. Patrz przykład 22 na str. 183$e$$e$

Zdarza się tak, że , gdzie jest zmienną losową, która jest zdefiniowana w następujący sposób. Jest to minimalna liczba taka, że i są liczbami losowymi z rozkładu równomiernego na . Piękne, prawda ?!ξ 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]$

Ponieważ jest to ćwiczenie, nie jestem pewien, czy fajnie jest dla mnie opublikować rozwiązanie (dowód) tutaj :) Jeśli chcesz sam to udowodnić, oto wskazówka: rozdział nazywa się „Chwile”, które powinny wskazywać jesteś we właściwym kierunku.

Jeśli chcesz wdrożyć go samodzielnie, nie czytaj dalej!

Jest to prosty algorytm do symulacji Monte Carlo. Narysuj jednolity losowy, a następnie kolejny i tak dalej, aż suma przekroczy 1. Liczba losowanych losów jest twoją pierwszą próbą. Powiedzmy, że masz:

 0.0180
 0.4596
 0.7920

Potem twoja pierwsza próba sprawiła, że ​​3. Kontynuuj te próby, a zauważysz, że średnio dostajesz .$e$

Poniżej znajduje się kod MATLAB, wynik symulacji i histogram.

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)

Wynik i 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

AKTUALIZACJA: Zaktualizowałem swój kod, aby pozbyć się szeregu wyników próbnych, aby nie zajmował pamięci RAM. Wydrukowałem również oszacowanie PMF.

Aktualizacja 2: Oto moje rozwiązanie Excel. Umieść przycisk w programie Excel i połącz go z następującym makrem 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

Wprowadź liczbę prób, na przykład 1000, w komórce D1 i kliknij przycisk. Oto jak powinien wyglądać ekran po pierwszym uruchomieniu:

AKTUALIZACJA 3: Silverfish zainspirował mnie do innej drogi, nie tak eleganckiej jak pierwsza, ale wciąż fajnej. Obliczył objętości n-simpleksów przy użyciu sekwencji Sobola .

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

Przypadkowo napisał pierwszą książkę o metodzie Monte Carlo, którą przeczytałem w szkole średniej. Moim zdaniem jest to najlepsze wprowadzenie do metody.

AKTUALIZACJA 4:

Silverfish w komentarzach sugerował prostą implementację formuły Excel. Taki wynik uzyskuje się po jego podejściu po około 1 milionie losowych liczb i 185 000 prób:

Oczywiście jest to znacznie wolniejsze niż wdrożenie Excel VBA. Zwłaszcza jeśli zmodyfikujesz mój kod VBA, aby nie aktualizować wartości komórek w pętli i zrobisz to dopiero po zebraniu wszystkich statystyk.

AKTUALIZACJA 5

Xi'an za rozwiązanie nr 3 jest ściśle powiązany (lub nawet taki sam w pewnym sensie jako komentarz na JWG w wątku). Trudno powiedzieć, kto wpadł na ten pomysł jako pierwszy Forsythe lub Gnedenko. Oryginalna edycja Gnedenko z 1950 roku w języku rosyjskim nie zawiera sekcji Problemy w rozdziałach. Tak więc nie mogłem znaleźć tego problemu na pierwszy rzut oka, gdzie jest w późniejszych wydaniach. Może został dodany później lub zakopany w tekście.

Jak skomentowałem w odpowiedzi Xi'ana, podejście Forsythe'a wiąże się z innym interesującym obszarem: rozkładem odległości między pikami (ekstrema) w losowych sekwencjach (IID). Średnia odległość zdarza się wynosić 3. Sekwencja dolna w podejściu Forsythe'a kończy się na dole, więc jeśli będziesz kontynuować próbkowanie, w pewnym momencie dostaniesz kolejne dno, a potem inne. Możesz śledzić odległość między nimi i zbudować rozkład.

Sugeruję, by głosować na odpowiedź Aksakala. Jest bezstronny i opiera się wyłącznie na metodzie generowania odchyleń jednostkowych jednostek.

Moja odpowiedź może być dowolnie sprecyzowana, ale nadal jest stronnicza od prawdziwej wartości .$e$

Odpowiedź Xi'ana jest poprawna, ale myślę, że jej zależność od funkcji lub sposobu generowania losowych odchyleń Poissona jest nieco okrągła, gdy celem jest przybliżenie .$\log$$e$

Oszacowanie przez ładowanie początkowe$e$

Zamiast tego rozważ procedurę ładowania. Jeden ma dużą liczbę obiektów które są rysowane z zamianą na wielkość próbki . Przy każdym losowaniu prawdopodobieństwo nie rysowania określonego obiektu wynosi i jest takich losowań. Prawdopodobieństwo pominięcia określonego obiektu we wszystkich losowaniach wynosin i 1 - n - 1 n p = ( 1 - $n$$n$$i$$1-n^{-1}$$n$$p=(1-\frac{1}{n})^n.$

Ponieważ zakładam, że wiemy, że

więc możemy też napisać

Oznacza to, że nasze oszacowanie można znaleźć poprzez oszacowanie prawdopodobieństwa pominięcia konkretnej obserwacji w bootstrapie replikuje w wielu takich replikacjach - tj. występowania obiektu w bootstrapie.m B j$p$$m$$B_j$$i$

Są dwa źródła błędów w tym przybliżeniu. Skończone zawsze oznacza, że ​​wyniki są przybliżone, tj. Oszacowanie jest stronnicze. Dodatkowo będzie oscylować wokół prawdziwej wartości, ponieważ jest to symulacja.$n$$\hat{p}$

Uważam to podejście za urocze, ponieważ student lub inna osoba z wystarczającą ilością rzeczy do zrobienia może przybliżyć za pomocą talii kart, stosu małych kamieni lub innych dostępnych przedmiotów, w tym samym stylu, co osoba może oszacować używając kompasu, prostej krawędzi i drobinek piasku. Myślę, że to fajne, kiedy matematykę można oddzielić od nowoczesnych udogodnień, takich jak komputery.$e$$\pi$

Wyniki

Przeprowadziłem kilka symulacji dla różnej liczby replik ładowania początkowego. Standardowe błędy są szacowane w normalnych odstępach czasu.

Zauważ, że wybór liczby obiektów ładowanych początkowo ustanawia absolutną górną granicę dokładności wyników, ponieważ procedura Monte Carlo szacuje a zależy tylko od . Ustawienie niepotrzebnie dużej wartości tylko obciąży komputer, albo dlatego, że potrzebujesz tylko „przybliżonego” przybliżenia do albo dlatego, że odchylenie zostanie zatłoczone przez wariancję z powodu Monte Carlo. Te wyniki są dla a jest dokładne z dokładnością do trzeciego miejsca po przecinku.p p n n e n = 10 3 p - 1 ≈ $n$$p$$p$$n$$n$$e$$n=10^3$$p^{-1}\approx e$

Ten wykres pokazuje, że wybór ma bezpośrednie i głębokie konsekwencje dla stabilności w . Niebieska linia przerywana pokazuje a czerwona linia pokazuje . Zgodnie z oczekiwaniami, zwiększenie wielkości próby daje coraz dokładniejsze szacunki . s t e $m$$\hat{p}$$p$$e$$\hat{p}$

Napisałem do tego żenująco długi skrypt R. Sugestie dotyczące ulepszeń można przesłać na odwrocie rachunku za 20 USD.

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')

Rozwiązanie 1:

Dla Poissona dystrybucja Dlatego, jeśli , co oznacza, że ​​możesz oszacować za pomocą symulacji Poissona. Symulacje Poissona można uzyskać z generatora rozkładu wykładniczego (jeśli nie w najbardziej efektywny sposób).P$\mathcal{P}(\lambda)$ X ∼ P ( 1 ) P ( X = 0 ) = P ( X = 1 ) = e - 1 e -

Uwaga 1: Jak omówiono w komentarzach, jest to dość skomplikowany argument, ponieważ symulowanie z rozkładu Poissona lub równoważnie rozkładu wykładniczego może być trudne do wyobrażenia bez udziału log lub funkcji exp ... Ale potem W. Huber doszedł do uratowanie tej odpowiedzi za pomocą najbardziej eleganckiego rozwiązania opartego na zamówionych mundurach. Co jest jednak przybliżeniem , ponieważ rozkład jednolitych odstępów jest Beta , co sugeruje, że który jest zbieżny do jako$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$rośnie do nieskończoności. Oprócz tego, że odpowiada na komentarze, generator wykładniczy von Neumanna z 1951 r. Używa tylko generacji jednorodnych.

Rozwiązanie 2:

Innym sposobem uzyskania reprezentacji stałej jako całki jest przypomnienie, że gdy to który jest również rozkładem . Dlatego Drugie podejście do przybliżenia przez Monte Carlo ma zatem symulować pary normalne i monitorować częstotliwość razy . W pewnym sensie jest to przeciwieństwo aproksymacji Monte Carlo związane z częstotliwością razy ...$e$

Rozwiązanie 3:

Mój kolega z Uniwersytetu w Warwick, M. Pollock, wskazał inne przybliżenie Monte Carlo zwane metodą Forsythe'a : chodzi o to, by uruchomić sekwencję ujednoliconych generacji aż do . Oczekiwanie na odpowiednią regułę zatrzymania, , która jest liczbą przypadków, w których jednorodna sekwencja spadła, wynosi wtedy podczas gdy prawdopodobieństwo, że jest nieparzyste, wynosi ! ( Metoda Forsythe'a faktycznie ma na celu symulację z dowolnej gęstości formy , stąd jest bardziej ogólna niż przybliżenie i .)u n + 1 > u n N e N e - 1 exp G ( x ) e $u_1,u_2,...$$u_{n+1}>u_{n}$$N$$e$$N$$e^{-1}$$\exp G(x)$$e$$e^{-1}$

Jest to dość równoległe do podejścia Gnedenko zastosowanego w odpowiedzi Aksakala , więc zastanawiam się, czy jedno można wyprowadzić z drugiego. Przynajmniej oba mają ten sam rozkład z masą prawdopodobieństwadla wartości .$1/n!$$n$

Szybka implementacja metody Forsythego R polega na rezygnacji z dokładnego wykonywania sekwencji mundurów na rzecz większych bloków, co pozwala na równoległe przetwarzanie:

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)

To nie rozwiązanie ... tylko szybki komentarz, który jest za długi na pole komentarza.

Aksakal

Aksakal opublikował rozwiązanie, w którym obliczamy oczekiwaną liczbę standardowych rysunków Uniform, które należy pobrać, tak aby ich suma przekroczyła 1. W Mathematica moim pierwszym sformułowaniem było:

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

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

EDYCJA: Właśnie się z tym pobawiłem, a następujący kod (ta sama metoda - także w Mma - tylko inny kod) jest około 10 razy szybszy:

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

Xian / Whuber

Whuber zasugerował szybki fajny kod do symulacji rozwiązania Xian 1:

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

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

który zauważa, że ​​jest 20 razy szybszy niż pierwszy kod (lub około dwa razy szybciej niż nowy kod powyżej).

Dla zabawy pomyślałem, że byłoby interesujące sprawdzić, czy oba podejścia są tak skuteczne (w sensie statystycznym). W tym celu wygenerowałem 2000 oszacowań e przy użyciu:

• Metoda Aksakala: dane A.
• Metoda Xiana 1 przy użyciu kodu Whuber: dataB

... oba w Mathematica . Poniższy diagram porównuje nieparametryczne oszacowanie gęstości jądra wynikowych zestawów danych dataA i dataB.

Tak więc, chociaż kod Whubera (czerwona krzywa) jest około dwa razy szybszy, metoda nie wydaje się być „niezawodna”.

Metoda wymagająca bezbożnej ilości próbek

$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}$

$\sqrt{2}$$\sqrt{2\pi}$

Metoda wymagająca bardzo niewielu próbek, ale powodująca bezbożną ilość błędów numerycznych

Całkowicie głupia, ale bardzo skuteczna odpowiedź na podstawie mojego komentarza:

$X \sim \text{uniform}(-1, 1)$$Y_n = | (\bar x)^n|$$\hat e = (1 - Y_n)^{-1/Y_n}$

To zbiegnie się bardzo szybko, ale również napotka ekstremalny błąd numeryczny.

$Y_n$$1/Y_n$$n \rightarrow \infty$$Y_n$$Y_n = 0$$e$

Oto inny sposób, w jaki można to zrobić, choć jest on dość wolny. Nie rości sobie pretensji do skuteczności, ale oferuję tę alternatywę w duchu kompletności.

$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$

Implementacja w R: Metodę można zaimplementować przy Rużyciu runifdo generowania jednolitych wartości. Kod jest następujący:

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$