Aproximado usando a simulação de Monte Carlo

35

Estive estudando a simulação de Monte Carlo recentemente e a tenho usado para aproximar constantes como (círculo dentro de um retângulo, área proporcional).π

No entanto, não consigo pensar em um método correspondente para aproximar o valor de e [número de Euler] usando a integração de Monte Carlo.

Você tem alguma indicação de como isso pode ser feito?

statisticsnewbie12345
fonte
7
Existem muitas, muitas, muitas maneiras de fazer isso. Que isso seja possível pode se tornar evidente ao contemplar o que o Rcomando 2 + mean(exp(-lgamma(ceiling(1/runif(1e5))-1)))faz. (Se o uso da função Gamma do log o incomoda, substitua-o por 2 + mean(1/factorial(ceiling(1/runif(1e5))-2)), que usa apenas adição, multiplicação, divisão e truncamento e ignore os avisos de estouro.) O que poderia ser de maior interesse seria simulações eficientes : você pode minimizar o número de etapas computacionais necessárias para estimar e com precisão?
whuber
4
Que pergunta deliciosa! Estou ansioso para ler as respostas dos outros. Uma maneira de você realmente chamar a atenção para essa pergunta - talvez outra meia dúzia de respostas - seria revisar a pergunta e pedir respostas eficientes , como sugere o whuber. Isso é como catnip para usuários de CV.
Sycorax diz Restabelecer Monica
11
@EngrStudent Não sei se o analógico geométrico existe para e . Simplesmente não é uma quantidade geométrica natural (com trocadilhos) como π .
Aksakal
6
@Aksakal e é uma quantidade excepcionalmente geométrica. No nível mais elementar, aparece naturalmente em expressões para áreas relacionadas à hipérbole. Em um nível um pouco mais avançado, está intimamente conectado a todas as funções periódicas, incluindo funções trigonométricas, cujo conteúdo geométrico é óbvio. O verdadeiro desafio aqui é que é tão fácil simular valores relacionados ao e !
whuber
2
@StatsStudent: e por si só não é interessante. No entanto, se isso levar a estimadores imparciais de quantidades como
exp{0xf(y)dG(y)}
isso poderá ser mais útil para os algoritmos MCMC.
Xian

Respostas:

34

A maneira simples e elegante de estimar por Monte Carlo é descrita neste artigo . O artigo é realmente sobre ensino . Portanto, a abordagem parece perfeitamente adequada ao seu objetivo. A ideia é baseada em um exercício de um popular livro russo sobre teoria das probabilidades de Gnedenko. Ver ex.22 na p.183eee

Isso acontece para que , onde é uma variável aleatória definida da seguinte maneira. É o número mínimo de tal que e são números aleatórios de distribuição uniforme em . Bonito, não é ?!ξ n n i = 1 r i > 1 r i [ 0 , 1 ]E[ξ]=eξni=1nri>1ri[0,1]

Como é um exercício, não tenho certeza se é legal postar a solução (prova) aqui :) Se você quiser provar você mesmo, aqui está uma dica: o capítulo é chamado de "Momentos", que deve apontar você na direção certa.

Se você deseja implementá-lo, não leia mais!

Este é um algoritmo simples para simulação de Monte Carlo. Desenhe um uniforme aleatório, depois outro e assim sucessivamente até que a soma exceda 1. O número de randoms sorteados é o seu primeiro teste. Digamos que você tenha:

 0.0180
 0.4596
 0.7920

Em seguida, seu primeiro teste foi renderizado 3. Continue fazendo esses testes e notará que, em média, você obtém .e

Seguem o código MATLAB, o resultado da simulação e o histograma.

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)

O resultado e o 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

insira a descrição da imagem aqui

ATUALIZAÇÃO: Atualizei meu código para livrar-se da matriz de resultados de testes, para que não ocupe RAM. Também imprimi a estimativa do PMF.

Atualização 2: Aqui está minha solução do Excel. Coloque um botão no Excel e vincule-o à seguinte macro 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

Digite o número de tentativas, como 1000, na célula D1 e clique no botão Aqui, como a tela deve ficar após a primeira execução:

insira a descrição da imagem aqui

ATUALIZAÇÃO 3: Silverfish me inspirou de outra maneira, não tão elegante quanto a primeira, mas ainda assim legal. Ele calculou os volumes de n-simplex usando sequências 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

Por coincidência, ele escreveu o primeiro livro sobre o método Monte Carlo que li no ensino médio. É a melhor introdução ao método na minha opinião.

ATUALIZAÇÃO 4:

O Silverfish nos comentários sugeriu uma implementação simples da fórmula do Excel. Esse é o tipo de resultado que você obtém com a abordagem dele após cerca de 1 milhão de números aleatórios e 185 mil tentativas:

insira a descrição da imagem aqui

Obviamente, isso é muito mais lento que a implementação do Excel VBA. Especialmente, se você modificar meu código VBA para não atualizar os valores das células dentro do loop, e só o fará quando todas as estatísticas forem coletadas.

ATUALIZAÇÃO 5

A solução 3 de Xi'an está intimamente relacionada (ou até a mesma em algum sentido, conforme o comentário de jwg no tópico). É difícil dizer quem teve a ideia primeiro em Forsythe ou Gnedenko. A edição original de Gnedenko em 1950, em russo, não possui seções de problemas nos capítulos. Portanto, à primeira vista, não foi possível encontrar esse problema em edições posteriores. Talvez tenha sido adicionado mais tarde ou enterrado no texto.

Como comentei na resposta de Xi'an, a abordagem de Forsythe está ligada a outra área interessante: a distribuição de distâncias entre picos (extremos) em seqüências aleatórias (IID). A distância média passa a ser 3. A sequência descendente na abordagem de Forsythe termina com um fundo, portanto, se você continuar amostrando, obterá outro fundo em algum momento, depois outro etc. Você pode rastrear a distância entre eles e criar a distribuição.

Aksakal
fonte
Uau isso é legal! Você poderia adicionar um parágrafo ou dois explicando por que isso funciona?
Sycorax diz Restabelecer Monica
7
(+1) Brilhante! A resposta merece a nota mais alta, pois depende apenas de simulações uniformes. E não usa nenhuma aproximação, exceto a de Monte Carlo. O fato de ele se conectar a Gnedenko é mais uma vantagem.
Xian
2
Legal! Aqui está Mathematica código para mesma, como um one-liner:
Mean[Table[ Length[NestWhileList[(Random[]+#) &, Random[], #<1&]], {10^6}]]
wolfies
4
@wolfies A seguinte tradução direta da Rsolução que postei na resposta de Xi'an é vinte vezes mais rápida:n=10^6; 1. / Mean[UnitStep[Differences[Sort[RandomReal[{0, n}, n + 1]]] - 1]]
whuber
11
Eu publiquei o "por que o significa ?" pergunta como uma pergunta em si mesma ; Eu suspeito que minha solução de esboço (que é o que imediatamente veio à mente como a visualização "óbvia" do problema) não é necessariamente a maneira que os estudantes russos pretendiam fazer isso! Portanto, soluções alternativas seriam muito bem-vindas. e
Silverfish
19

Sugiro revogar a resposta de Aksakal. É imparcial e depende apenas de um método para gerar desvios uniformes da unidade.

Minha resposta pode ser arbitrariamente precisa, mas ainda assim é desviada do verdadeiro valor de .e

A resposta de Xi'an está correta, mas acho que sua dependência da função ou de uma maneira de gerar desvios aleatórios de Poisson é um pouco circular quando o objetivo é aproximar .eloge

Estimando pelo Bootstrappinge

Em vez disso, considere o procedimento de inicialização. Um deles tem um grande número de objetos que são desenhados com substituição para um tamanho de amostra de . A cada sorteio, a probabilidade de não desenhar um objeto específico é , e existem tais sorteios. A probabilidade de um objeto específico ser omitido de todos os desenhos én i 1 - n - 1 n p = ( 1 - 1nni1n1np=(11n)n.

Porque eu suponho que sabemos que

exp(1)=limn(11n)n

para que também possamos escrever

exp(1)p^=i=1mIiBjm

Ou seja, nossa estimativa de é encontrada estimando-se a probabilidade de uma observação específica ser omitida de bootstrap replica em muitas dessas repetições - ou seja, a fração de ocorrências do objeto nas bootstraps.m B j ipmBji

Existem duas fontes de erro nessa aproximação. finito sempre significa que os resultados são aproximados, ou seja, a estimativa é enviesada. Além disso, flutuará em torno do valor verdadeiro, porque esta é uma simulação.pnp^

Acho essa abordagem um tanto encantadora, porque um estudante de graduação ou outra pessoa com pouco o que fazer poderia se aproximar usando um baralho de cartas, uma pilha de pedras pequenas ou quaisquer outros itens disponíveis, na mesma linha que uma pessoa poderia estimar usando uma bússola, uma régua e alguns grãos de areia. Eu acho que é legal quando a matemática pode ser separada das conveniências modernas, como computadores.πeπ

Resultados

Realizei várias simulações para vários números de replicações de inicialização. Os erros padrão são estimados usando intervalos normais.

Note-se que a escolha de o número de objetos que estão sendo bootstrapped fixa um máximo absoluto sobre a precisão dos resultados porque o procedimento Monte Carlo está estimando e depende apenas . Definir como desnecessariamente grande apenas sobrecarregará seu computador, porque você só precisa de uma aproximação "aproximada" de ou porque o viés será inundado pela variação devido ao Monte Carlo. Estes resultados são para e são precisas para a terceira casa decimal.p p n n e n = 10 3 p - 1enppnnen=103p-1 1e

Este gráfico mostra que a escolha de tem conseqüências diretas e profundas para a estabilidade em . A linha tracejada azul mostra e a linha vermelha mostra . Como esperado, aumentar o tamanho da amostra produz estimativas cada vez mais precisas . p p e pmp^pep^insira a descrição da imagem aqui

Eu escrevi um roteiro R embaraçosamente longo para isso. Sugestões para melhoria podem ser enviadas no verso de uma nota de US $ 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')
Sycorax diz restabelecer Monica
fonte
4
+1 Faz muito sentido. Alguma chance de compartilhar seu código se você o escreveu?
Antoni Parellada
2
Mesmo que isso pode ser arbitrariamente precisas, em última análise, não é satisfatória porque ele só simula uma aproximação para , em vez de em si. eee
whuber
11
Certo. Você terminaria com uma chamada replicada dentro de outra, que é essencialmente a mesma que temos agora.
Sycorax diz Restabelecer Monica
11
@whuber Eu realmente não vejo a distinção entre uma aproximação arbitrariamente precisa de uma aproximação arbitrariamente precisa de , e uma aproximação arbitrariamente precisa de . eee
GTC
11
@jwg Além de ser conceitualmente importante, também é praticamente importante porque implementar uma aproximação a uma aproximação requer acompanhar o quão precisa é cada uma das duas aproximações. Mas eu teria que concordar que, quando ambas as aproximações são aceitáveis, a abordagem geral é boa.
whuber
14

Solução 1:

Para uma distribuição Poisson , Portanto, se , que significa que você pode estimar por uma simulação de Poisson. E as simulações de Poisson podem ser derivadas de um gerador de distribuição exponencial (se não da maneira mais eficiente).P ( X = k ) = λ kP(λ) X P ( 1 ) P ( X = 0 ) = P ( X = 1 ) = e - 1 e - 1

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

Observação 1: Como discutido nos comentários, esse é um argumento bastante complicado, pois pode ser difícil imaginar uma simulação de uma distribuição Poisson ou equivalentemente uma distribuição exponencial sem envolver um log ou uma função exp ... Mas, então, W. Huber chegou ao resgate desta resposta com uma solução mais elegante baseada em uniformes encomendados. Porém, o que é uma aproximação , pois a distribuição de um espaçamento uniforme é um Beta , implicando que que converge para comoU(i:n)U(i1:n)B(1,n)

P(n{U(i:n)U(i1:n)}1)=(11n)n
e1ncresce até o infinito. Como outro aparte que responde aos comentários, o gerador exponencial de von Neumann de 1951 usa apenas gerações uniformes.

Solução 2:

Outra maneira de obter uma representação da constante como uma integral é lembrar que, quando depois que também é uma distribuição . Portanto, Uma segunda abordagem para aproximar por Assim, Monte Carlo simula pares normais e monitora a frequência dos tempos . Em certo sentido, é o oposto da aproximação de Monte Carlo de relacionada à frequência dos tempos ...X 1 , X 2 iid ~ N ( 0 , 1 ) ( X 2 1 + X 2 2 ) ~ χ - 1 e ( X 1 ,e

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

Solução 3:

M. Pollock, meu colega da Universidade Warwick, apontou outra aproximação de Monte Carlo chamada método de Forsythe : a idéia é executar uma sequência de gerações uniformes até . A expectativa da regra de parada correspondente, , que é o número de vezes que a sequência uniforme foi desativada é então enquanto a probabilidade de ser ímpar é ! ( O método de Forsythe na verdade visa simular a partir de qualquer densidade da forma , portanto, é mais geral do que aproximar-se de e .)você1 1,você2,...un+1>unNeNe1expG(x)ee1

Isso é bastante paralelo à abordagem de Gnedenko usada na resposta de Aksakal , então me pergunto se uma pode ser derivada da outra. No mínimo, ambos têm a mesma distribuição com massa de probabilidadepara o valor .n1/n!n

Uma rápida implementação em R do método de Forsythe é deixar de seguir com precisão a sequência de uniformes em favor de blocos maiores, o que permite o processamento 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)
Xi'an
fonte
12
Desde que se saiba como fazer a simulação de Poisson sem saber . e
Glen_b -Reinstala Monica
5
Se eu ligar para o gerador R rpoiss (), posso fingir que não conheço . Mais seriamente, você pode gerar variáveis ​​exponenciais [usando uma função vez de ] até que a soma exceda e o número resultante menos um seja Poisson . E ( 1 ) log e 1 P ( 1 )eE(1 1)registroe1 1P(1 1)
Xian
5
Computar equivale a comput , uma vez que são inversas. Você pode evitar o cálculo de qualquer função de várias maneiras. Aqui está uma solução baseada na sua primeira resposta: Ela usa apenas aritmética elementar. expregistroexpn <- 1e5; 1/mean(n*diff(sort(runif(n+1))) > 1)
whuber
3
Eu acredito que o método de Forsythe é o mesmo que o de Gnedenko. Escolher um uniforme, de modo que seja menor que 1, é o mesmo que escolher menor que , e se tivermos êxito, é distribuído uniformemente condicionalmente entre e 0.n x i x n 1 - n - 1 x i 1 - n x i 1 - n - 1 x ixnnxixn1n1xi1nxi1n1xEu
jwg 12/02/16
3
Eu não estava ciente da abordagem de Forsythe. No entanto, está ligado a algo muito interessante. Se em vez de parar em você continuar a amostragem, em seguida, a expectativa da distância de para o próximo final é exatamente 3.nn+1 1n
Aksakal
7

Não é uma solução ... apenas um comentário rápido que é muito longo para a caixa de comentários.

Aksakal

Aksakal postou uma solução em que calculamos o número esperado de desenhos uniformes padrão que devem ser feitos, de modo que sua soma exceda 1. No Mathematica , minha primeira formulação foi:

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

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

EDIT: Apenas joguei rápido com isso, e o código a seguir (mesmo método - também no Mma - apenas código diferente) é cerca de 10 vezes mais rápido:

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

Xian / Whuber

Whuber sugeriu código legal rápido para simular a solução 1 de Xian:

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

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

que ele observa é 20 vezes mais rápido que o primeiro código (ou cerca de duas vezes mais rápido que o novo código acima).

Apenas por diversão, pensei que seria interessante ver se as duas abordagens são tão eficientes (no sentido estatístico). Para fazer isso, gerei 2000 estimativas de e usando:

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

... ambos no Mathematica . O diagrama a seguir contrasta uma estimativa não-paramétrica da densidade do kernel dos conjuntos de dados dataA e dataB resultantes.

insira a descrição da imagem aqui

Portanto, embora o código do whuber (curva vermelha) seja duas vezes mais rápido, o método não parece ser 'confiável'.

wolfies
fonte
Uma linha vertical no local do valor real melhoraria bastante essa imagem.
Sycorax diz Restabelecer Monica
11
É uma observação muito interessante, obrigado. Como a meia largura será dimensionada quadraticamente com o tamanho da simulação e a meia largura do método de Xi'an é aproximadamente o dobro do método de Aksakal, executar quatro vezes mais iterações as tornará igualmente precisas. A questão de quanto esforço é necessário em cada iteração permanece: se uma iteração do método de Xi'an exigir menos de um quarto do esforço, esse método ainda será mais eficiente.
whuber
11
Acredito que a situação se torna clara quando você compara o número de realizações de variáveis ​​aleatórias necessárias nos dois métodos, em vez do valor nominal de . n
whuber
11
@whuber escreveu: running four times as many iterations will make them equally accurate///// ..... Acabei de jogar rápido com isso: aumentando o número de pontos de amostra usados ​​no método 1 de Xian de para 6 x (ou seja, 6 vezes o número de pontos) produz uma curva semelhante à de Aksaksal. 10 6106106
wolfies
11
Bem feito com o código - será difícil melhorar muito nisso.
whuber
2

Método que requer uma quantidade ímpia de amostras

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

N(0 0,1 1)ex

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

22π

Método que requer muito poucas amostras, mas causa uma quantidade injusta de erro numérico

Uma resposta completamente boba, mas muito eficiente, com base em um comentário que fiz:

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

Isso convergirá muito rápido, mas também ocorrerá um erro numérico extremo.

Yn1 1/YnnYnYn=0 0e

Cliff AB
fonte
2
e
11
@ whuber: eu não usei o Box-Muller devido à transformação de log necessária muito diretamente para exponencial no meu livro. Eu teria permitido reflexivamente cos e pecado, mas isso foi apenas porque eu havia esquecido uma análise complexa por um momento, um ponto tão bom.
Cliff AB
11
n1 1n2ϕ(2)n1 1n2en2n1 1
2

Aqui está outra maneira de fazer isso, embora seja bastante lento. Não reivindico eficiência, mas ofereço essa alternativa no espírito de perfeição.

nvocê1 1,,vocênIID U(0 0,1 1)e

E(Eu(vocêEu1 1/e)vocêEu)=1 1/e1 1dvocêvocê=1

evocê(1 1)você(n)

Sn(k)1 1nEu=1 1k1 1você(Eu)para todos k=1 1,..,n.

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

e^2você(m)+você(m+1 1).

1 1/eee

Implementação em R: O método pode ser implementado no Ruso runifpara gerar valores uniformes. O código é o seguinte:

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

Restabelecer Monica
fonte