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 [número de Euler] usando a integração de Monte Carlo.
Você tem alguma indicação de como isso pode ser feito?
simulation
monte-carlo
algorithms
random-generation
numerical-integration
statisticsnewbie12345
fonte
fonte
R
comando2 + mean(exp(-lgamma(ceiling(1/runif(1e5))-1)))
faz. (Se o uso da função Gamma do log o incomoda, substitua-o por2 + 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 estimarRespostas:
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.183ee e
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 ξ n ∑ni=1ri>1 ri [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:
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.
O resultado e o histograma:
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:
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:
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 .
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:
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.
fonte
Mean[Table[ Length[NestWhileList[(Random[]+#) &, Random[], #<1&]], {10^6}]]
R
soluçã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]]
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 .elog e
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 - 1n n i 1−n−1 n p=(1−1n)n.
Porque eu suponho que sabemos que
para que também possamos escrever
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 ip m Bj i
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.pn p^
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 - 1 ≈ en p p n n e n = 103 p- 1≈ e
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 pm p^ p e p^
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.
fonte
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
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
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, u2, . . . vocên+1>un N e N e−1 expG(x) e e−1
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:
fonte
n <- 1e5; 1/mean(n*diff(sort(runif(n+1))) > 1)
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:
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:
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:
... 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.
Portanto, embora o código do whuber (curva vermelha) seja duas vezes mais rápido, o método não parece ser 'confiável'.
fonte
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 6Método que requer uma quantidade ímpia de amostras
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:
Isso convergirá muito rápido, mas também ocorrerá um erro numérico extremo.
fonte
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.
Implementação em R: O método pode ser implementado no
R
usorunif
para gerar valores uniformes. O código é o seguinte:fonte