Aproximado de $e$ utilizando la Simulación de Monte Carlo

Question

He estado mirando la simulación de Monte Carlo recientemente, y han estado utilizando a la aproximación de las constantes tales como $\pi$ (círculo dentro de un rectángulo, proporcionadas área).

Sin embargo, soy incapaz de pensar en un método para aproximar el valor de $e$ utilizando la simulación de Monte Carlo.

Tienes alguna sugerencia sobre cómo puede hacerse esto?

Answer 1

La sencilla y elegante forma de estimar la $e$ por Monte Carlo se describe en este documento. El papel es en realidad acerca de la enseñanza de $e$. Por lo tanto, el enfoque parece perfectamente adecuado para su objetivo. La idea se basa en un ejercicio de una popular ruso libro de texto sobre la teoría de la probabilidad por Gnedenko. Ver ex.22 en la p.183

Sucede así que el $E[\xi]=e$ donde $\xi$ es una variable aleatoria que se define como sigue. Es el número mínimo de $n$ tal que $\sum_{i=1}^n r_i>1$ $r_i$ son números aleatorios de una distribución uniforme en $[0,1]$. Hermoso, ¿no?!

Ya que es un ejercicio, no estoy seguro de si es bueno para mí para publicar la solución (prueba) aquí :) Si quieres probar a sí mismo, he aquí un consejo: el capítulo se llama "Momentos", que debería apuntar en la dirección correcta.

Si desea aplicar a ti mismo, entonces no hay que leer más!

Este es un sencillo algoritmo para la simulación de Monte Carlo. Dibujar un aleatorio uniforme, luego otro y así sucesivamente hasta que la suma sea superior a 1. El número de randoms dibujado es su primer juicio. Digamos que tengo:

 0.0180
 0.4596
 0.7920

Luego de su primer juicio prestados 3. Seguir haciendo estas pruebas, y te darás cuenta de que, en promedio, usted consigue $e$.

Código de MATLAB, resultados de simulación y el histograma de seguir.

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)

El resultado y el histograma:

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

ACTUALIZACIÓN: He actualizado mi código para deshacerse de la matriz de resultados de los ensayos para que no tome RAM. Yo también impreso el PMF estimación.

Actualización 2: Aquí está mi solución Excel. Poner un botón en Excel y enlace a la siguiente macro de 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

Introduzca el número de ensayos, tales como 1000, en la celda D1, y haga clic en el botón. He aquí cómo la pantalla se verá como después de la primera carrera:

ACTUALIZACIÓN 3: El pececillo de plata inspirado a mí de otra manera, no tan elegante como el primero pero todavía fresco. Se calculan los volúmenes de n-simplexes utilizando Sobol secuencias.

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

Casualmente él escribió el primer libro sobre el método de Monte Carlo me vuelva a leer en la escuela secundaria. Es la mejor introducción al método, en mi opinión.

ACTUALIZACIÓN 4:

El pececillo de plata en los comentarios sugirieron una simple fórmula de Excel implementación. Este es el tipo de resultado que se obtiene con su estrategia después de que sobre el total de 1 millón de números aleatorios y 185K ensayos:

Obviamente, esto es mucho más lento que el VBA de Excel implementación. Especialmente, si usted modificar mi código de VBA para no actualizar los valores de las celdas dentro del bucle, y sólo lo hacen una vez que todas las estadísticas son recogidos.

Answer 2

Sugiero upvoting Aksakal la respuesta. Es imparcial y se basa sólo en un método de generación de la unidad de uniforme se desvía.

Mi respuesta puede ser hecho de forma arbitraria precisa, pero todavía está sesgada lejos del verdadero valor de $e$.

Xi'an, la respuesta es correcta, pero creo que su dependencia en el $\log$ o de la función de una forma de generación aleatoria de Poisson se desvía un poco de la circular cuando el propósito es aproximar $e$.

La estimación de $e$ por el Arranque

En su lugar, considere el procedimiento de bootstrapping. Uno tiene un gran número de objetos de $n$ que se dibujan con la sustitución de un tamaño de muestra de $n$. En cada sorteo, la probabilidad de no dibujar un objeto en particular $i$$1-n^{-1}$, y $n$ estos sorteos. La probabilidad de que un objeto en particular, se omite de todos los sorteos es $p=(1-\frac{1}{n})^n.$

Porque estoy suponiendo que sabemos que $$\exp(-1)=\lim_{n\to\infty}\left(1-\frac{1}{n}\right)^n$$

por lo tanto, también puede escribir $$\exp(-1)\approx \hat{p}=\sum_{i=1}^m\frac{\mathbb{I}_{i\in B_j}}{m}$$

That is, our estimate of $p$ is found by estimating the probability that a specific observation is omitted from $m$ bootstrap replicates $B_j$ across many such replicates -- i.e. the fraction of occurrences of object $i$ in the bootstraps.

There are two sources of error in this approximation. Finite $n$ will always mean that the results are approximate, i.e. the estimate is biased. Additionally, $\hat{p}$ will fluctuate around the true value because this is a simulation.

I find this approach somewhat charming because an undergraduate or another person with sufficiently little to do could approximate $e$ using a deck of cards, a pile of small stones, or any other items at hand, in the same vein as a person could estimate $\pi$ using a compass, a straight-edge and some grains of sand. I think it's neat when mathematics can be divorced from modern conveniences like computers.

Results

I conducted several simulations for various number of bootstrap replications. Standard errors are estimated using normal intervals.

Note that the choice of $n$ the number of objects being bootstrapped sets an absolute upper limit on the accuracy of the results because the Monte Carlo procedure is estimating $p$ and $p$ depends only on $n$. Setting $n$ to be unnecessarily large will just encumber your computer, either because you only need a "rough" approximation to $e$ or because the bias will be swamped by variance due to the Monte Carlo. These results are for $n=10^3$ and $p^{-1}\aprox e$ is accurate to the third decimal.

This plot shows that the choice of $m$ has direct and profound consequences for the stability in $\hat{p}$. The blue dashed line shows $p$ and the red line shows $e$. As expected, increasing the sample size produces ever-more accurate estimates $\hat{p}$.

Escribí una embarazosa largo R script para esta. Sugerencias para la mejora puede ser presentado en el reverso de un billete de $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')

Answer 3

Solución 1:

Para una distribución de Poisson $\mathcal{P}(\lambda)$ distribución, $$\mathbb{P}(X=k)=\frac{\lambda^k}{k!}\,e^{-\lambda}$$Therefore, if $X\sim\mathcal{P}(1)$, $$\mathbb{P}(X=0)=\mathbb{P}(X=1)=e^{-1}$$which means you can estimate $e^{-1}$ Poisson de simulación.

Como se discute en los comentarios, es más bien un enrevesado argumento desde la simulación de una distribución de Poisson o, equivalentemente, una Exponencial la distribución es difícil de imaginar sin la participación de un registro o una exp la función... Pero entonces W. Huber vino con una solución más elegante basado en ordenadas los uniformes. Que es una aproximación sin embargo, dado que la distribución de un uniforme espaciado $U_{(i:n)}-U_{(i-1:n)}$ es una Beta $\mathfrak{B}(1,n)$, lo que implica que $$\mathbb{P}(n\{U_{(i:n)}-U_{(i-1:n)}\}\ge 1)=\left(1-\frac{1}{n}\right)^n$$which converges to $e^{-1}$ as $$ n crece hasta el infinito.

Solución 2:

Otra forma de lograr una representación de la constante de $e$ como una integral es recordar que, cuando $$X_1,X_2\stackrel{\text{iid}}{\sim}\mathfrak{N}(0,1)$$ then $$(X_1^2+X_2^2)\sim\chi^2_1$$ which is also an $\mathcal{E}(1/2)$ de distribución. Por lo tanto, $$\mathbb{P}(X_1^2+X_2^2\ge 2)=1-\{1-\exp(-2/2)\}=e^{-1}$$ Una segunda aproximación a la aproximación de $e$ por Monte Carlo es, pues, para simular las condiciones normales de los pares de $(X_1,X_2)$ y monitorear la frecuencia de los tiempos de $X_1^2+X_2^2\ge 2$. En un sentido es el opuesto de la Monte Carlo aproximación de $\pi$ con relación a la frecuencia de los tiempos de $X_1^2+X_2^2<1$...

Answer 4

No es una solución ... sólo un breve comentario de que es demasiado largo para el cuadro de comentario.

Aksakal

Aksakal publicado una solución en la que se calcula el número esperado de los Uniformes de los dibujos que se deben tomar, tales que la suma exceda de 1. En Mathematica, mi primera formulación fue:

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

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

EDIT: acaban de tener un juego rápido con esto, y el siguiente código (el mismo método también en las Mma, simplemente es diferente de código) es aproximadamente 10 veces más rápido:

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

Xian / Whuber

Whuber ha sugerido frío rápido código para simular Xian solución 1:

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

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

lo que se observa es 20 veces más rápido que el primer código (o aproximadamente dos veces tan rápido como el nuevo código de arriba).

Sólo por diversión, pensé que sería interesante ver si ambos enfoques son tan eficientes (en el sentido estadístico). Para ello, se generó el año 2000 las estimaciones de correo utilizando:

• Aksakal del método: dataA
• Xian en el método 1 uso de whuber código: datosb

... tanto en Mathematica. En el siguiente diagrama se contrasta un test no paramétrico de estimación de densidad de kernel de la resultante de dataA y datosb conjuntos de datos.

Así, mientras que whuber del código (curva roja) es el doble de rápido, el método no parece ser tan "fiable".

Answer 5

Método que requiere un impío cantidad de muestras

En primer lugar usted necesita para ser capaz de muestra de una distribución normal. Suponiendo que usted va a excluir el uso de la función de $f(x) = e^x$, o buscar mesas derivada de esa función, se puede producir aproximado de muestras de la distribución normal a través de la CLT. Por ejemplo, si usted puede degustar de un uniforme(0,1) la distribución,$\frac{ \bar x \sqrt{12}}{ \sqrt{n}} \dot \sim N(0,1)$. Como se ha señalado por whuber, para tener la última estimación enfoque de $e$ como el tamaño de la muestra enfoques $\infty$, sería necesario que el número de uniforme de las muestras que se utilicen enfoques $\infty$ como el tamaño de la muestra enfoques infinito.

Ahora, si se puede degustar de una distribución normal, con grandes muestras suficientes, usted puede obtener estimaciones consistentes de la densidad de $N(0,1)$. Esto se puede hacer con los histogramas o núcleo suavizadores (pero tenga cuidado de no utilizar un núcleo Gaussiano para seguir su no$e^{x}$!!). Para conseguir sus estimaciones de densidad para ser coherente, usted tendrá que dejar su df (número de contenedores en el histograma, la inversa de la ventana para hacer más suaves) enfoque infinito, pero más lento que el tamaño de la muestra.

Así que ahora, con una gran cantidad de potencia de cálculo, puede calcular la densidad de una $N(0,1)$, es decir,$\hat \phi(x)$. Desde $\phi(\sqrt(2) ) = (2 \pi)^{-1/2} e^{-1}$, la estimación de $e = \hat \phi(\sqrt{2}) \sqrt{2 \pi}$.

Si quieres ir totalmente frutos secos, incluso se puede estimar el $\sqrt{2}$ $\sqrt{2\pi}$ utilizando los métodos descritos anteriormente.

Método que requiere muy pocas muestras, pero causando un impío cantidad de error numérico

Un completo tonto, pero muy eficiente, responder a partir de un comentario que hice:

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

Esto converge muy rápido, pero también se ejecutan en la más extrema de error numérico.

whuber señaló que este usa la función de potencia, lo que normalmente se llama a la función exp. Esto podría ser evadido por la discretización $Y_n$, de tal manera que $1/Y_n$ es un número entero y el poder podría ser reemplazado con la multiplicación repetida. Sería necesario que como $n \rightarrow \infty$, el de la discretización de $Y_n$ obtener más fino y más fino,y la discretización tendría que excluir $Y_n = 0$. Con esto, el estimador de la teoría (es decir, el mundo en el que el error numérico no existe) convergería a $e$, y bastante rápido!