Ungefähre mit Hilfe von Monte - Carlo - Simulation

35

Ich habe mir kürzlich die Monte-Carlo-Simulation angesehen und sie verwendet, um Konstanten wie (Kreis in einem Rechteck, proportionale Fläche) anzunähern.π

Ich kann mir jedoch keine entsprechende Methode vorstellen, um den Wert von [Eulers Zahl] mithilfe der Monte-Carlo-Integration zu approximieren .e

Haben Sie Hinweise, wie dies getan werden kann?

statisticnewbie12345
quelle
7
Es gibt viele, viele, viele Möglichkeiten, dies zu tun. Dass dies so ist, könnte sich zeigen, wenn man darüber nachdenkt, was der RBefehl 2 + mean(exp(-lgamma(ceiling(1/runif(1e5))-1)))tut. (Wenn Sie die Gamma-Protokollfunktion stören, ersetzen Sie diese durch 2 + mean(1/factorial(ceiling(1/runif(1e5))-2))Addition, Multiplikation, Division und Trunkierung und ignorieren Sie die Überlaufwarnungen.) Was von größerem Interesse sein könnte, wären effiziente Simulationen: Können Sie die Anzahl minimieren ? Rechenschritte erforderlich, um mit einer bestimmten Genauigkeit abzuschätzen ? e
Whuber
4
Was für eine entzückende Frage! Ich freue mich darauf, die Antworten anderer zu lesen. Eine Möglichkeit, die Aufmerksamkeit auf diese Frage zu lenken - vielleicht ein weiteres halbes Dutzend Antworten - wäre, die Frage zu überarbeiten und nach effizienten Antworten zu fragen , wie Whuber vorschlägt. Das ist wie eine Katzenminze für CV-Benutzer.
Sycorax sagt Reinstate Monica
1
@EngrStudent Ich bin mir nicht sicher, ob das geometrische Analog für existiert . Es ist einfach keine natürliche (Wortspiel beabsichtigte) geometrische Größe wie . πeπ
Aksakal
6
@Aksakal ist eine außergewöhnlich geometrische Größe. Auf der elementarsten Ebene erscheint es natürlich in Ausdrücken für Bereiche, die mit Hyperbeln zusammenhängen. Auf einer etwas fortgeschritteneren Ebene ist es eng mit allen periodischen Funktionen verbunden, einschließlich der trigonometrischen Funktionen, deren geometrischer Inhalt offensichtlich ist. Die eigentliche Herausforderung dabei ist, dass es so einfach ist , bezogene Werte zu simulieren ! ee
Whuber
2
@StatsStudent: an sich ist nicht interessant. Wenn dies jedoch zu unverzerrten Schätzern von Größen wie kann sich dies für MCMC-Algorithmen als am nützlichsten erweisen. exp { x 0 f ( y ) d G ( y ) }e
exp{0xf(y)dG(y)}
Xi'an

Antworten:

34

Die einfache und elegante Art, nach Monte Carlo zu schätzen, wird in diesem Artikel beschrieben . In der Arbeit geht es eigentlich um das Unterrichten von . Daher scheint der Ansatz perfekt zu Ihrem Ziel zu passen. Die Idee basiert auf einer Übung aus einem beliebten russischen Lehrbuch zur Wahrscheinlichkeitstheorie von Gnedenko. Siehe Beispiel 22 auf S.183eee

Es passiert also, dass , wobei eine Zufallsvariable ist, die wie folgt definiert ist. Es ist die minimale Anzahl von so dass und Zufallszahlen aus der Gleichverteilung von . Schön, nicht wahr ?!ξ n n i = 1 r i > 1 r i [ 0 , 1 ]E[ξ]=eξni=1nri>1ri[0,1]

Da es sich um eine Übung handelt, bin ich mir nicht sicher, ob es für mich cool ist, die Lösung (Beweis) hier zu posten :) Wenn Sie es selbst beweisen möchten, hier ein Tipp: Das Kapitel heißt "Momente", was darauf hindeuten sollte Sie in die richtige Richtung.

Wenn Sie es selbst implementieren möchten, lesen Sie nicht weiter!

Dies ist ein einfacher Algorithmus für die Monte-Carlo-Simulation. Ziehe einen einheitlichen Zufallsgenerator, dann einen weiteren und so weiter, bis die Summe 1 übersteigt. Die Anzahl der gezogenen Zufälle ist dein erster Versuch. Angenommen, Sie haben:

 0.0180
 0.4596
 0.7920

Dann wurde Ihre erste Testversion gerendert. 3. Führen Sie diese Testversionen weiter aus, und Sie werden feststellen, dass Sie im Durchschnitt .e

MATLAB-Code, Simulationsergebnis und das Histogramm folgen.

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)

Das Ergebnis und das Histogramm:

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

Bildbeschreibung hier eingeben

UPDATE: Ich habe meinen Code aktualisiert, um die Reihe der Testergebnisse zu entfernen, sodass kein RAM benötigt wird. Ich habe auch die PMF-Schätzung ausgedruckt.

Update 2: Hier ist meine Excel-Lösung. Fügen Sie eine Schaltfläche in Excel ein und verknüpfen Sie sie mit dem folgenden VBA-Makro:

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

Geben Sie die Anzahl der Versuche (z. B. 1000) in die Zelle D1 ein und klicken Sie auf die Schaltfläche. So sollte der Bildschirm nach dem ersten Lauf aussehen:

Bildbeschreibung hier eingeben

UPDATE 3: Silverfish hat mich auf eine andere Art inspiriert, nicht so elegant wie die erste, aber trotzdem cool. Es berechnete das Volumen von n-Simplexen unter Verwendung von Sobol- Sequenzen.

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

Zufällig schrieb er das erste Buch über die Monte-Carlo-Methode, das ich in der High School gelesen hatte. Meiner Meinung nach ist dies die beste Einführung in die Methode.

UPDATE 4:

Silverfish schlug in Kommentaren eine einfache Implementierung von Excel-Formeln vor. Dies ist die Art von Ergebnis, die Sie mit seiner Herangehensweise nach insgesamt 1 Million Zufallszahlen und 185K-Versuchen erhalten:

Bildbeschreibung hier eingeben

Dies ist offensichtlich viel langsamer als die Implementierung von Excel VBA. Insbesondere, wenn Sie meinen VBA-Code so ändern, dass die Zellenwerte in der Schleife nicht aktualisiert werden und erst dann aktualisiert werden, wenn alle Statistiken erfasst wurden.

UPDATE 5

Xi'ans Lösung Nr. 3 ist eng miteinander verwandt (oder in gewisser Weise sogar mit jwgs Kommentar im Thread identisch). Es ist schwer zu sagen, wer zuerst Forsythe oder Gnedenko auf die Idee kam. Gnedenkos Originalausgabe von 1950 in russischer Sprache enthält keine Abschnitte mit Problemen in Kapiteln. So konnte ich dieses Problem auf den ersten Blick nicht finden, wo es in späteren Ausgaben ist. Vielleicht wurde es später hinzugefügt oder im Text vergraben.

Wie ich in Xi'ans Antwort bemerkte, ist Forsythes Ansatz mit einem anderen interessanten Bereich verbunden: der Verteilung der Abstände zwischen Spitzen (Extrema) in zufälligen (IID) Sequenzen. Der mittlere Abstand ist zufällig 3. Die Abwärtssequenz in Forsythes Ansatz endet mit einem Boden. Wenn Sie also mit dem Abtasten fortfahren, erhalten Sie irgendwann einen anderen Boden, dann einen anderen usw. Sie könnten den Abstand zwischen ihnen verfolgen und die Verteilung aufbauen.

Aksakal
quelle
Wow, das ist cool! Könnten Sie einen oder zwei Absätze hinzufügen, in denen erläutert wird, warum dies funktioniert?
Sycorax sagt Reinstate Monica
7
(+1) Genial! Die Antwort verdient die Bestnote, da sie sich nur auf einheitliche Simulationen stützt. Und verwendet keine andere Annäherung als die von Monte Carlo. Dass es sich wieder mit Gnedenko verbindet, ist ein weiterer Vorteil.
Xi'an
2
Cool! Hier ist der Mathematica- Code für dasselbe als
Mean[Table[ Length[NestWhileList[(Random[]+#) &, Random[], #<1&]], {10^6}]]
Einzeiler
4
@wolfies Die folgende direkte Übersetzung der RLösung, die ich in Xi'ans Antwort gepostet habe, ist zwanzigmal schneller:n=10^6; 1. / Mean[UnitStep[Differences[Sort[RandomReal[{0, n}, n + 1]]] - 1]]
whuber
1
Ich habe die "Warum ist der Mittelwert ?" Frage als eine Frage für sich ; Ich vermute, dass meine Skizzenlösung (was mir sofort als "offensichtliche" Visualisierung des Problems in den Sinn kam) nicht unbedingt so ist, wie es die russischen Studenten beabsichtigten! Daher wären alternative Lösungen sehr willkommen. e
Silverfish
19

Ich schlage vor, Aksakals Antwort zu unterstützen. Es ist unvoreingenommen und stützt sich nur auf ein Verfahren zur Erzeugung einheitlicher Abweichungen.

Meine Antwort kann willkürlich präzisiert werden, ist aber immer noch vom wahren Wert von .e

Xi'ans Antwort ist richtig, aber ich denke, dass seine Abhängigkeit von der Funktion oder einer Methode zur Erzeugung von Poisson-Zufallsabweichungen ein wenig kreisförmig ist, wenn der Zweck darin besteht, e anzunähern .loge

Schätzung von durch Bootstrappinge

Betrachten Sie stattdessen das Bootstrapping-Verfahren. Man hat eine große Anzahl von Objekten die beim Ersetzen auf eine Stichprobengröße von n gezeichnet werden . Bei jeder Ziehung beträgt die Wahrscheinlichkeit , ein bestimmtes Objekt i nicht zu zeichnen, 1 - n - 1 , und es gibt n solche Ziehungen. Die Wahrscheinlichkeit, dass ein bestimmtes Objekt in allen Ziehungen weggelassen wird, ist p = ( 1 - 1nni1n1np=(11n)n.

exp(1)=limn(11n)n

so können wir auch schreiben

exp(1)p^=i=1mIiBjm

Das heißt, unsere Schätzung von ergibt sich aus der Schätzung der Wahrscheinlichkeit, dass eine bestimmte Beobachtung in Bootstrap-Replikaten über viele solcher Replikate hinweg unterbleibt - dh dem Anteil des Auftretens von Objekt in den Bootstraps.m B j ipmBji

Bei dieser Annäherung gibt es zwei Fehlerquellen. Endliches immer, dass die Ergebnisse ungefähr sind, dh die Schätzung ist verzerrt. Außerdem schwankt um den wahren Wert, da dies eine Simulation ist.pnp^

Ich finde diesen Ansatz etwas charmant , weil ein Student oder eine andere Person mit ausreichend wenig zu tun annähern könnte mit einem Deck von Karten, einen Haufen von kleinen Steinen oder andere Gegenständen in der Hand, in der gleichen Richtung wie eine Person schätzen , könnte mit einem Kompass, einer geraden Kante und ein paar Sandkörnern. Ich finde es gut, wenn Mathematik von modernen Annehmlichkeiten wie Computern getrennt werden kann.πeπ

Ergebnisse

Ich habe mehrere Simulationen für verschiedene Anzahl von Bootstrap-Replikationen durchgeführt. Standardfehler werden in normalen Intervallen geschätzt.

Beachten Sie, dass die Auswahl von die Anzahl der Objekte, die gebootet werden, eine absolute Obergrenze für die Genauigkeit der Ergebnisse darstellt, da die Monte-Carlo-Prozedur schätzt und nur von abhängt . Wenn Sie unnötig groß einstellen, belastet dies nur Ihren Computer, entweder weil Sie nur eine "grobe" Annäherung an benötigen oder weil die Abweichung aufgrund des Monte-Carlo-Effekts überlastet wird. Diese Ergebnisse sind für und auf die dritte Dezimalstelle genau.p p n n e n = 10 3 p - 1enppnnen=103p1e

Diese Darstellung zeigt, dass die Wahl von direkte und tiefgreifende Konsequenzen für die Stabilität in . Die blaue gestrichelte Linie zeigt und die rote Linie zeigt . Wie erwartet, erzeugt die Probengröße zu erhöhen immer genauere Schätzungen . p p e pmp^pep^Bildbeschreibung hier eingeben

Ich habe dafür ein peinlich langes R-Skript geschrieben. Verbesserungsvorschläge können auf der Rückseite einer 20-Dollar-Rechnung eingereicht werden.

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')
Sycorax sagt Reinstate Monica
quelle
4
+1 Es macht sehr viel Sinn. Gibt es eine Chance, dass Sie Ihren Code teilen können, wenn Sie ihn geschrieben haben?
Antoni Parellada
2
ee
1
Sicher. Sie würden einfach mit einem Wiederholungsaufruf in einem anderen enden, der im Wesentlichen derselbe ist, den wir jetzt haben.
Sycorax sagt Reinstate Monica
1
ee
1
@jwg Das ist nicht nur konzeptionell wichtig, sondern auch praktisch wichtig, da für die Implementierung einer Annäherung an eine Annäherung protokolliert werden muss, wie genau die beiden Annäherungen sind. Aber ich muss zustimmen, dass, wenn beide Annäherungen akzeptabel sind, der Gesamtansatz in der Tat in Ordnung ist.
Whuber
14

Lösung 1:

P(λ)

P(X=k)=λkk!eλ
XP(1)
P(X=0)=P(X=1)=e1
e1

U(i:n)U(i1:n)B(1,n)

P(n{U(i:n)U(i1:n)}1)=(11n)n
e1n

Lösung 2:

e

X1,X2iidN(0,1)
(X12+X22)χ12
E(1/2)
P(X12+X222)=1{1exp(2/2)}=e1
e(X1,X2)X12+X222πX12+X22<1

Lösung 3:

u1,u2,...un+1>unNeNe1expG(x)ee1

1/n!n

Eine schnelle R-Implementierung der Forsythe-Methode besteht darin, die Reihenfolge der Uniformen zugunsten größerer Blöcke genau einzuhalten, was eine parallele Verarbeitung ermöglicht:

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)
Xi'an
quelle
12
e
5
eE(1)loge1P(1)
5
logexpn <- 1e5; 1/mean(n*diff(sort(runif(n+1))) > 1)
3
xnnxixn1n1xi1nxi1n1xi
3
n+1n
7

Keine Lösung ... nur ein kurzer Kommentar, der für das Kommentarfeld zu lang ist.

Aksakal

Aksakal veröffentlichte eine Lösung, in der wir die erwartete Anzahl der zu erstellenden einheitlichen Standardzeichnungen so berechnen, dass ihre Summe 1 übersteigt. In Mathematica lautete meine erste Formulierung:

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

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

EDIT: Hatte gerade ein schnelles Spiel damit und der folgende Code (gleiche Methode - auch in Mma - nur anderer Code) ist ungefähr 10-mal schneller:

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

Xian / Whuber

Whuber hat schnellen coolen Code vorgeschlagen, um Xians Lösung 1 zu simulieren:

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

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

was er feststellt, ist 20 mal schneller als der erste Code (oder etwa doppelt so schnell wie der neue Code oben).

Nur zum Spaß dachte ich, es wäre interessant zu sehen, ob beide Ansätze (im statistischen Sinne) so effizient sind. Zu diesem Zweck habe ich 2000 Schätzungen für e erstellt, in denen Folgendes verwendet wurde:

  • Aksakals Methode: dataA
  • Xians Methode 1 unter Verwendung des Whuber-Codes: dataB

... beides in Mathematica . Das folgende Diagramm stellt eine nichtparametrische Schätzung der Kerneldichte der resultierenden Datensätze dataA und dataB gegenüber.

Bildbeschreibung hier eingeben

Obwohl Whubers Code (rote Kurve) etwa doppelt so schnell ist, scheint die Methode nicht so zuverlässig zu sein.

Wolfies
quelle
Eine vertikale Linie am Ort des wahren Wertes würde dieses Bild erheblich verbessern.
Sycorax sagt Reinstate Monica
1
Es ist eine sehr interessante Beobachtung, danke. Da die Halbwertsbreite quadratisch mit der Größe der Simulation skaliert und die Halbwertsbreite von Xi'ans Methode etwa doppelt so groß ist wie die von Aksakal, werden sie durch Ausführen von viermal so vielen Iterationen gleich genau. Die Frage, wie viel Aufwand für jede Iteration erforderlich ist, bleibt offen: Wenn eine Iteration der Xi'an-Methode weniger als ein Viertel des Aufwands in Anspruch nimmt, ist diese Methode immer noch effizienter.
Whuber
1
n
1
running four times as many iterations will make them equally accurate106106
1
Gut gemacht mit dem Code - es wird schwierig sein, das zu verbessern.
Whuber
2

Methode, die eine gottlose Menge an Proben erfordert

f(x)=exx¯12n˙N(0,1)e

N(0,1)ex

N(0,1)ϕ^(x)ϕ((2))=(2π)1/2e1e=ϕ^(2)2π

Wenn Sie total verrückt werden wollen, können Sie sogar schätzen22π

Methode, die nur sehr wenige Stichproben erfordert, aber eine gottlose Menge an numerischen Fehlern verursacht

Eine völlig dumme, aber sehr effiziente Antwort, die auf einem Kommentar basiert, den ich abgegeben habe:

Xuniform(1,1)Yn=|(x¯)n|e^=(1Yn)1/Yn

Dies wird sehr schnell konvergieren , aber auch auf extreme numerische Fehler stoßen.

Yn1/YnnYnYn=0e

Cliff AB
quelle
2
e
1
@whuber: Ich habe den Box-Muller aufgrund der erforderlichen Log-Transformation nicht zu direkt in exponentiell in meinem Buch verwendet. Ich hätte Lattich und Sünde reflexartig zugelassen, aber das war nur, weil ich für einen Moment die komplexe Analyse vergessen hatte, so ein guter Punkt.
Cliff AB
1
n1n2ϕ(2)n1n2en2n1
2

Hier ist eine andere Möglichkeit, es zu tun, obwohl es ziemlich langsam ist. Ich mache keinen Anspruch auf Effizienz, biete diese Alternative aber im Sinne der Vollständigkeit an.

nU1,,UnIID U(0,1)e

E(I(Ui1/e)Ui)=1/e1duu=1.

eu(1)u(n)

Sn(k)1ni=1k1u(i)for all k=1,..,n.

mmin{k|S(k)1}1/ee

e^2u(m)+u(m+1).

1/eee

Implementierung in R: Die Methode kann implementiert werden, Rindem runifeinheitliche Werte generiert werden. Der Code lautet wie folgt:

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

e

Setzen Sie Monica wieder ein
quelle