He estado mirando la simulación de Monte Carlo recientemente, y la he estado usando para aproximar constantes como (círculo dentro de un rectángulo, área proporcional).$\pi$

Sin embargo, no puedo pensar en un método correspondiente para aproximar el valor de [número de Euler] utilizando la integración de Monte Carlo.$e$

¿Tiene alguna sugerencia sobre cómo se puede hacer esto?

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

Sucede de modo que , donde es una variable aleatoria que se define de la siguiente manera. Es el número mínimo de tal que y son números aleatorios de distribución uniforme en . Hermoso, ¿no?ξ 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]$

Dado que es un ejercicio, no estoy seguro de si es bueno para mí publicar la solución (prueba) aquí :) Si desea probarlo usted mismo, aquí hay un consejo: el capítulo se llama "Momentos", que debería señalar usted en la dirección correcta

Si desea implementarlo usted mismo, ¡no siga leyendo!

Este es un algoritmo simple para la simulación de Monte Carlo. Dibuja un uniforme al azar, luego otro y así sucesivamente hasta que la suma exceda 1. El número de randoms extraídos es tu primer intento. Digamos que tienes:

 0.0180
 0.4596
 0.7920

Luego, su primera prueba se procesó. 3. Siga haciendo estas pruebas y notará que, en promedio, obtiene .$e$

El código MATLAB, el resultado de la simulación y el histograma siguen.

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: Actualicé mi código para deshacerme de la matriz de resultados de prueba para que no tome RAM. También imprimí la estimación de PMF.

Actualización 2: Aquí está mi solución de Excel. Ponga un botón en Excel y vincúlelo 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

Ingrese el número de pruebas, como 1000, en la celda D1 y haga clic en el botón. Así es como debería verse la pantalla después de la primera ejecución:

ACTUALIZACIÓN 3: Silverfish me inspiró de otra manera, no tan elegante como la primera pero aún genial. Calculó los volúmenes de n-simplexes usando secuencias de 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

Casualmente escribió el primer libro sobre el método de Monte Carlo que leí en la escuela secundaria. Es la mejor introducción al método en mi opinión.

ACTUALIZACIÓN 4:

Silverfish en los comentarios sugirió una implementación simple de la fórmula de Excel. Este es el tipo de resultado que obtienes con su enfoque después de aproximadamente 1 millón de números aleatorios y 185K pruebas:

Obviamente, esto es mucho más lento que la implementación de Excel VBA. Especialmente, si modifica mi código VBA para no actualizar los valores de las celdas dentro del bucle, y solo lo hace una vez que se recopilan todas las estadísticas.

ACTUALIZACIÓN 5

De Xi'an solución # 3 está estrechamente relacionado (o incluso el mismo en cierto sentido como un comentario de por jwg en el hilo). Es difícil decir quién se le ocurrió la idea primera Forsythe o Gnedenko. Edición original de 1950 Gnedenko en ruso no tiene secciones problemas en los capítulos. Por lo tanto, no pude encontrar este problema a primera vista, donde se encuentra en ediciones posteriores. Tal vez fue agregado más tarde o enterrado en el texto.

Como comenté en la respuesta de Xi'an, el enfoque de Forsythe está vinculado a otra área interesante: la distribución de distancias entre picos (extremos) en secuencias aleatorias (IID). La distancia media es 3. La secuencia descendente en el enfoque de Forsythe termina con un fondo, por lo que si continúa muestreando obtendrá otro fondo en algún momento, luego otro, etc. Puede rastrear la distancia entre ellos y construir la distribución.

Sugiero votar la respuesta de Aksakal. Es imparcial y se basa únicamente en un método para generar desviaciones uniformes de la unidad.

Mi respuesta puede hacerse arbitrariamente precisa, pero todavía está sesgada lejos del verdadero valor de .$e$

La respuesta de Xi'an es correcta, pero creo que su dependencia de la función o de una forma de generar desviaciones aleatorias de Poisson es un poco circular cuando el propósito es aproximar .$\log$$e$

La estimación de por Bootstrapping$e$

Por el contrario, considerar el procedimiento de bootstrapping. Uno tiene un gran número de objetos que se dibujan con el reemplazo a un tamaño de muestra de . En cada sorteo, la probabilidad de no dibujar un objeto particular es , y hay tales empates. La probabilidad de que un objeto en particular se omite de todos los sorteos esn i 1 - n - 1 n p = ( 1 - $n$$n$$i$$1-n^{-1}$$n$$p=(1-\frac{1}{n})^n.$

Porque supongo que sabemos que

entonces también podemos escribir

Es decir, nuestra estimación de se encuentra al estimar la probabilidad de que se omita una observación específica de bootstrap replica en muchas de tales réplicas, es decir, la fracción de ocurrencias del objeto en los bootstraps.$p$$m$$B_j$$i$

Hay dos fuentes de error en esta aproximación. Finito siempre significará que los resultados son aproximados, es decir, la estimación está sesgada. Además, fluctuará alrededor del valor verdadero porque es una simulación.$n$$\hat{p}$

Encuentro este enfoque algo encantador porque un estudiante universitario u otra persona con poco que hacer podría aproximarse a usando una baraja de cartas, un montón de piedras pequeñas o cualquier otro elemento a la mano, en la misma línea que una persona podría estimar usando una brújula, un borde recto y algunos granos de arena. Creo que es bueno cuando las matemáticas se pueden divorciar de las comodidades modernas como las computadoras.$e$$\pi$

Resultados

Realicé varias simulaciones para varias réplicas de bootstrap. Los errores estándar se estiman utilizando intervalos normales.

Tenga en cuenta que la elección de el número de objetos que se están bootstrapped conjuntos de un límite superior absoluto en la precisión de los resultados debido a que el procedimiento de Monte Carlo es la estimación de y sólo depende de . Establecer para que sea innecesariamente grande solo afectará a su computadora, ya sea porque solo necesita una aproximación "aproximada" de o porque el sesgo se verá afectado por la variación debido al Monte Carlo. Estos resultados son para y es exacto al tercer decimal.$n$$p$$p$$n$$n$$e$$n=10^3$$p^{-1}\approx e$

Este gráfico muestra que la elección de tiene consecuencias directas y profundas para la estabilidad en . La línea discontinua azul muestra y la línea roja muestra . Como se esperaba, aumentar el tamaño de la muestra produce estimaciones cada vez más precisas . $m$$\hat{p}$$p$$e$$\hat{p}$

Escribí un guión R vergonzosamente largo para esto. Se pueden enviar sugerencias de mejora al reverso de una factura 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')

Solución 1:

Para una distribución de Poisson , Por lo tanto, si , que significa que puede estimar por una simulación de Poisson. Y las simulaciones de Poisson se pueden derivar de un generador de distribución exponencial (si no de la manera más eficiente).$\mathcal{P}(\lambda)$

Observación 1: Como se discutió en los comentarios, este es un argumento bastante complicado ya que simular a partir de una distribución de Poisson o, de manera equivalente, una distribución exponencial puede ser difícil de imaginar sin involucrar un registro o una función exp ... Pero entonces W. Huber llegó a la rescate de esta respuesta con la solución más elegante basada en uniformes ordenados. Sin embargo, es una aproximación , ya que la distribución de un espaciado uniforme es un Beta , lo que implica que que converge a como$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$crece hasta el infinito Como otro lado que responde a los comentarios, el generador exponencial de 1951 de von Neumann solo usa generaciones uniformes.

Solución 2:

Otra forma de conseguir una representación de la constante como una integral es recordar que, cuando entonces que es también un distribución. Por lo tanto, Un segundo enfoque para aproximar por Monte Carlo es, pues, para simular pares normales y controlar la frecuencia de veces que . En cierto sentido, es lo opuesto a la aproximación de Monte Carlo relacionada con la frecuencia de veces ...$e$

Solución 3:

Mi colega de la Universidad de Warwick, M. Pollock, señaló otra aproximación de Montecarlo llamada método de Forsythe : la idea es ejecutar una secuencia de generaciones uniformes hasta . ¡La expectativa de la regla de detención correspondiente, , que es el número de veces que la secuencia uniforme bajó es entonces mientras que la probabilidad de que sea ​​impar es ! ( El método de Forsythe en realidad apunta a simular desde cualquier densidad de la forma , por lo tanto, es más general que aproximar y .)$u_1,u_2,...$$u_{n+1}>u_{n}$$N$$e$$N$$e^{-1}$$\exp G(x)$$e$$e^{-1}$

Esto es bastante paralelo al enfoque de Gnedenko utilizado en la respuesta de Aksakal , por lo que me pregunto si uno puede derivarse del otro. ¡Como mínimo, ambos tienen la misma distribución con probabilidad de masapor valor .$1/n!$$n$

Una implementación rápida de R del método de Forsythe es renunciar a seguir con precisión la secuencia de uniformes a favor de bloques más grandes, lo que permite el procesamiento paralelo:

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)

No es una solución ... solo un comentario rápido que es demasiado largo para el cuadro de comentarios.

Aksakal

Aksakal publicó una solución en la que calculamos el número esperado de dibujos uniformes estándar que deben tomarse, de modo que su suma exceda 1. En Mathematica , mi primera formulación fue:

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

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

EDITAR: Acabo de jugar rápidamente con esto, y el siguiente código (el mismo método, también en Mma, solo un código diferente) 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 un código rápido y rápido para simular la solución 1 de Xian:

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

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

que él nota es 20 veces más rápido que el primer código (o aproximadamente el doble de rápido que el nuevo código anterior).

Solo por diversión, pensé que sería interesante ver si ambos enfoques son tan eficientes (en un sentido estadístico). Para hacerlo, generé 2000 estimaciones de e usando:

• Método de Aksakal: dataA
• Método 1 de Xian usando el código whuber: dataB

... ambos en Mathematica . El siguiente diagrama contrasta una estimación de densidad de kernel no paramétrica de los conjuntos de datos dataA y dataB resultantes.

Entonces, mientras que el código de Whuber (curva roja) es aproximadamente el doble de rápido, el método no parece ser tan 'confiable'.

Método que requiere una cantidad impía de muestras

Primero debe poder tomar muestras de una distribución normal. Suponiendo que va a excluir el uso de la función , o buscar tablas derivadas de esa función, puede producir muestras aproximadas de la distribución normal a través del CLT. Por ejemplo, si puede tomar muestras de una distribución uniforme (0,1), entonces . Como señaló Whuber, para tener el enfoque de estimación final medida que el tamaño de la muestra se aproxima a , se requeriría que el número de muestras uniformes utilizadas se acerque a medida que el tamaño de la muestra se aproxima al infinito.$f(x) = e^x$$\frac{ \bar x \sqrt{12}}{ \sqrt{n}} \dot \sim N(0,1)$$e$$\infty$$\infty$

Ahora, si puede tomar muestras de una distribución normal, con muestras suficientemente grandes, puede obtener estimaciones consistentes de la densidad de . Esto se puede hacer con histogramas o suavizadores de kernel (¡pero tenga cuidado de no usar un kernel gaussiano para seguir su regla no !). Para que sus estimaciones de densidad sean consistentes, deberá dejar que su df (número de bins en el histograma, inverso de la ventana para suavizar) se acerque al infinito, pero más lento que el tamaño de la muestra.$N(0,1)$$e^{x}$

Entonces, ahora, con mucha potencia computacional, puede aproximar la densidad de un , es decir, . Como , su estimación para .$N(0,1)$$\hat \phi(x)$$\phi(\sqrt(2) ) = (2 \pi)^{-1/2} e^{-1}$$e = \hat \phi(\sqrt{2}) \sqrt{2 \pi}$

Si quiere volverse completamente loco, incluso puede estimar y usando los métodos que discutió anteriormente.$\sqrt{2}$$\sqrt{2\pi}$

Método que requiere muy pocas muestras, pero que causa una cantidad impía de error numérico

Una respuesta completamente tonta, pero muy eficiente, basada en un comentario que hice:

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

Esto convergerá muy rápido, pero también se encontrará con un error numérico extremo.

Whuber señaló que esto utiliza la función de potencia, que normalmente llama a la función exp. Esto podría discretizando , de modo que sea ​​un número entero y la potencia se pueda reemplazar con la multiplicación repetida. Sería necesario que como , la discretización de hiciera cada vez más fina, y la discretización tendría que excluir . Con esto, el estimador teóricamente (es decir, el mundo en el que no existe un error numérico) convergería en , ¡y bastante rápido!$Y_n$$1/Y_n$$n \rightarrow \infty$$Y_n$$Y_n = 0$$e$

Aquí hay otra forma de hacerlo, aunque es bastante lento. No pretendo ser eficiente, pero ofrezco esta alternativa en un espíritu de integridad.

Contra la respuesta de Xi'an a los fines de esta pregunta que usted puede generar y usar una secuencia de variables pseudoaleatorias uniformes y luego necesita estimar por algún método usando operaciones aritméticas básicas (es decir, no puede usar funciones logarítmicas o exponenciales ni ninguna distribución que use estas funciones). El presente método está motivado por un resultado simple que involucra variables aleatorias uniformes:$n$$U_1, \cdots , U_n \sim \text{IID U}(0,1)$$e$$^\dagger$

Estimación de usando este resultado:$e$ Primero ordenamos los valores de la muestra en orden descendente para obtener las estadísticas de orden y luego definimos las sumas parciales:$u_{(1)} \geqslant \cdots \geqslant u_{(n)}$

Ahora, dejemos que y luego estimar por interpolación de las variables uniformes ordenadas. Esto proporciona un estimador para dado por:$m \equiv \min \{ k | S(k) \geqslant 1 \}$$1/e$$e$

Este método tiene un ligero sesgo (debido a la interpolación lineal del punto de corte para ) pero es un estimador consistente para . El método puede implementarse con bastante facilidad, pero requiere la clasificación de valores, que es más computacionalmente intensivo que el cálculo determinista de . Este método es lento, ya que implica la clasificación de valores.$1/e$$e$$e$

Implementación en R: El método se puede implementar al Rusar runifpara generar valores uniformes. El código es el siguiente:

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]); }

La implementación de este código da convergencia al verdadero valor de , pero es muy lento en comparación con los métodos deterministas.$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$ Creo que queremos evitar cualquier método que haga uso de cualquier transformación que implique un exponencial o un logaritmo. Si podemos usar densidades que usan exponenciales en su definición, entonces es posible derivar de ellas algebraicamente usando una llamada de densidad.$e$