Je me suis intéressé récemment à la simulation de Monte Carlo et je l’utilise pour approcher des constantes telles que (cercle à l’intérieur d’un rectangle, zone proportionnelle).
Cependant, je suis incapable de penser à une méthode correspondante pour approximer la valeur de [nombre d'Euler] en utilisant l'intégration de Monte Carlo.
Avez-vous des indications sur la façon dont cela peut être fait?
simulation
monte-carlo
algorithms
random-generation
numerical-integration
statisticnewbie12345
la source
la source
R
commandement2 + mean(exp(-lgamma(ceiling(1/runif(1e5))-1)))
fait. (Si vous utilisez la fonction de consignation Gamma vous dérange, remplacez-la par2 + mean(1/factorial(ceiling(1/runif(1e5))-2))
, qui utilise uniquement l’ajout, la multiplication, la division et la troncature, et ignorez les avertissements de débordement.) Ce qui pourrait être plus intéressant serait des simulations efficaces : pouvez-vous réduire le nombre de étapes de calcul nécessaires pour estimerRéponses:
La méthode simple et élégante pour estimere par Monte Carlo est décrite dans cet article . Le papier est en fait sur l'enseignement e . Par conséquent, l'approche semble parfaitement adaptée à votre objectif. L'idée est basée sur un exercice d'un manuel populaire russe sur la théorie des probabilités de Gnedenko. Voir ex.22 p.183
Il se trouve queE[ξ]=e , où ξ est une variable aléatoire qui est définie comme suit. C'est le nombre minimum de n tel que ∑ni=1ri>1 et ri sont des nombres aléatoires de distribution uniforme sur [0,1] . Beau, n'est ce pas?
Puisqu'il s'agit d'un exercice, je ne suis pas sûr que ce soit cool pour moi de poster la solution (preuve) ici :) Si vous souhaitez le prouver vous-même, voici un conseil: le chapitre s'appelle "Moments", ce qui devrait indiquer vous dans la bonne direction.
Si vous voulez l'implémenter vous-même, alors ne lisez pas plus loin!
C'est un algorithme simple pour la simulation de Monte Carlo. Dessinez un tirage au sort uniforme, puis un autre et ainsi de suite jusqu'à ce que la somme dépasse 1. Le nombre de tirages au sort correspond à votre premier essai. Disons que vous avez:
Ensuite, votre premier essai rendu 3. Continuez à faire ces essais, et vous remarquerez qu'en moyenne, vous obteneze .
Le code MATLAB, le résultat de la simulation et l'histogramme suivent.
Le résultat et l'histogramme:
MISE À JOUR: J'ai mis à jour mon code pour supprimer le tableau des résultats des essais afin qu'il ne prenne pas de mémoire vive. J'ai également imprimé l'estimation du PMF.
Mise à jour 2: Voici ma solution Excel. Placez un bouton dans Excel et associez-le à la macro VBA suivante:
Entrez le nombre d'essais, par exemple 1000, dans la cellule D1 et cliquez sur le bouton. Voici à quoi l'écran devrait ressembler après la première utilisation:
UPDATE 3: Silverfish m'a inspiré d'une autre manière, pas aussi élégante que la première mais toujours cool. Il a calculé les volumes de n-simplexes à l'aide de séquences Sobol .
Par coïncidence, il a écrit le premier livre sur la méthode de Monte Carlo que j'ai lu au lycée. C'est la meilleure introduction à la méthode à mon avis.
MISE À JOUR 4:
Silverfish dans les commentaires a suggéré une implémentation simple de la formule Excel. C'est le genre de résultat obtenu avec son approche après environ 1 million de nombres aléatoires et 185 000 essais:
Évidemment, cela est beaucoup plus lent que l'implémentation Excel VBA. En particulier, si vous modifiez mon code VBA pour ne pas mettre à jour les valeurs de cellule dans la boucle et ne le faites que lorsque toutes les statistiques sont collectées.
MISE À JOUR 5
La solution n ° 3 de Xi'an est étroitement liée (ou même identique dans un sens, comme le dit le commentaire de Jwg dans le fil de discussion). Il est difficile de dire qui a eu cette idée en premier Forsythe ou Gnedenko. L'édition originale de 1950 de Gnedenko en russe ne comporte pas de sections Problèmes dans les chapitres. Donc, je ne pouvais pas trouver ce problème au premier abord où il se trouve dans les éditions ultérieures. Peut-être que cela a été ajouté plus tard ou enterré dans le texte.
Comme je l'ai commenté dans la réponse de Xi'an, l'approche de Forsythe est liée à un autre domaine intéressant: la distribution des distances entre les pics (extrema) dans des séquences aléatoires (IID). La distance moyenne dans l’approche de Forsythe se termine par un fond, donc si vous continuez à échantillonner, vous obtiendrez un autre fond à un moment donné, puis un autre, etc. Vous pouvez suivre la distance qui les sépare et construire la distribution.
la source
Mean[Table[ Length[NestWhileList[(Random[]+#) &, Random[], #<1&]], {10^6}]]
R
solution que j'ai publiée dans la réponse de Xi'an est vingt fois plus rapide:n=10^6; 1. / Mean[UnitStep[Differences[Sort[RandomReal[{0, n}, n + 1]]] - 1]]
Je suggère de voter pour la réponse d'Aksakal. Il est non biaisé et repose uniquement sur une méthode permettant de générer des écarts uniformes.
Ma réponse peut être arbitrairement précise, mais reste biaisée par rapport à la vraie valeur dee .
La réponse de Xi'an est correcte, mais je pense que sa dépendance à la fonction de ou à un moyen de générer des déviations aléatoires de Poisson est un peu circulaire lorsque le but est d'approximer e .log e
Estimation de par Bootstrappinge
Au lieu de cela, considérons la procédure d'amorçage. On a un grand nombre d'objets qui sont dessinés avec remplacement à une taille d'échantillon de n . A chaque tirage, la probabilité de ne pas dessiner un objet particulier i est 1 - n - 1 , et il y a n tirages. La probabilité qu’un objet particulier soit omis de tous les tirages est de p = ( 1 - 1n n i 1−n−1 n p=(1−1n)n.
Parce que je suppose que nous savons que
donc on peut aussi écrire
C'est-à-dire que notre estimation de est obtenue en estimant la probabilité qu'une observation spécifique soit omise dans m répliques bootstrap B j sur plusieurs répliques de ce type - c'est-à-dire la fraction d'occurrences de l'objet i dans les bootstraps.p m Bj i
Il y a deux sources d'erreur dans cette approximation. Fini signifie toujours que les résultats sont approximatifs, c'est-à-dire que l'estimation est biaisée. En outre, pn p^ fluctuera autour de la valeur réelle parce que c'est une simulation.
Je trouve cette approche un peu de charme , car un premier cycle ou d'une autre personne avec suffisamment peu pour ne pourraient approcher en utilisant un jeu de cartes, un tas de petites pierres ou d'autres éléments à portée de main, dans la même veine que une personne pourrait estimer π en utilisant une boussole, une règle et des grains de sable. Je pense que c'est bien quand les mathématiques peuvent être dissociées des commodités modernes comme les ordinateurs.e π
Résultats
J'ai effectué plusieurs simulations pour différents nombres de réplications bootstrap. Les erreurs-types sont estimées à l'aide d'intervalles normaux.
Notez que le choix de nombre d'objets à initialiser définit une limite supérieure absolue pour la précision des résultats, car la procédure de Monte Carlo estime que p et p ne dépendent que de n . Définir n comme étant inutilement grand encombrera simplement votre ordinateur, soit parce que vous n’avez besoin que d’une approximation "approximative" de e, soit parce que le biais sera submergé par la variance due au Monte Carlo. Ces résultats sont pour n = 10 3 et p - 1 ≈ e est précis à la troisième décimale.n p p n n e n=103 p−1≈e
Ce graphique montre que le choix des a des conséquences directes et profondes pour la stabilité en p . La ligne pointillée bleue indique p et la ligne rouge indique e . Comme prévu, l' augmentation de la taille de l' échantillon produit des estimations toujours plus précises p .m p^ p e p^
J'ai écrit un script R d'une longueur embarrassante pour cela. Les suggestions d'amélioration peuvent être soumises au verso d'un billet de 20 $.
la source
Solution 1:
Pour une distribution de Poisson , P ( X = k ) = λ kP( λ ) Donc, si X ∼ P ( 1 ) ,
P ( X = 0 ) = P ( X = 1 ) = e - 1, ce qui signifie que vous pouvez estimer e - 1 par une simulation de Poisson. Et les simulations de Poisson peuvent être dérivées d'un générateur de distribution exponentielle (si ce n'est de la manière la plus efficace).
Solution 2:
Solution 3:
A quick R implementation of Forsythe's method is to forgo following precisely the sequence of uniforms in favour of larger blocks, which allows for parallel processing:
la source
n <- 1e5; 1/mean(n*diff(sort(runif(n+1))) > 1)
It uses only elementary arithmetic.Pas une solution ... juste un commentaire rapide qui est trop long pour la zone de commentaire.
Aksakal
Aksakal a publié une solution dans laquelle nous calculons le nombre attendu de dessins uniformes standard à prendre, de sorte que leur somme dépasse 1. Dans Mathematica , ma première formulation était la suivante:
EDIT: Je viens d’avoir un jeu rapide avec cela, et le code suivant (même méthode - également en Mma - code tout différent) est environ 10 fois plus rapide:
Xian / Whuber
Whuber a suggéré un code rapide pour simuler la solution de Xian 1:
Version R:
n <- 1e5; 1/mean(n*diff(sort(runif(n+1))) > 1)
Version Mma:
n=10^6; 1. / Mean[UnitStep[Differences[Sort[RandomReal[{0, n}, n + 1]]] - 1]]
note qu'il est 20 fois plus rapide que le premier code (ou environ deux fois plus vite que le nouveau code ci-dessus).
Juste pour le plaisir, j'ai pensé qu'il serait intéressant de voir si les deux approches sont aussi efficaces (au sens statistique du terme). Pour ce faire, j'ai généré 2000 estimations de e en utilisant:
... both in Mathematica. The following diagram contrasts a nonparametric kernel density estimate of the resulting dataA and dataB data sets.
So, while whuber's code (red curve) is about twice as fast, the method does not appear to be as 'reliable'.
la source
running four times as many iterations will make them equally accurate
///// ..... Je viens d'avoir un jeu rapide avec ceci: augmenter le nombre de points d'échantillon utilisés dans la méthode 1 de Xian à partir deMéthode nécessitant une quantité impie d'échantillons
Vous devez d'abord pouvoir échantillonner à partir d'une distribution normale. En supposant que vous allez exclure l'utilisation de la fonctionF( x ) = eX ou recherchez des tables dérivées de cette fonction, vous pouvez produire des échantillons approximatifs de la distribution normale via le CLT. Par exemple, si vous pouvez échantillonner à partir d’une distribution uniforme (0,1), alorsX¯12√n√~˙N( 0 , 1 ) . Comme l'a souligné Whuber, adopter l'approche de l'estimation finalee à mesure que la taille de l'échantillon approche ∞ , il serait nécessaire que le nombre d'échantillons uniformes utilisés approche ∞ lorsque la taille de l'échantillon approche l'infini.
Maintenant, si vous pouvez échantillonner à partir d'une distribution normale, avec des échantillons suffisamment grands, vous pouvez obtenir des estimations cohérentes de la densité deN( 0 , 1 ) . Cela peut être fait avec des histogrammes ou des lisseurs de noyau (mais veillez à ne pas utiliser de noyau gaussien pour suivre votre non.eX règle!). Pour que vos estimations de densité soient cohérentes, vous devez laisser votre df (nombre de cases dans l'histogramme, inverse de la fenêtre pour obtenir une surface plus lisse) approcher l'infini, mais plus lentement que la taille de l'échantillon.
Alors maintenant, avec beaucoup de puissance de calcul, vous pouvez approximer la densité d'unN( 0 , 1 ) , c'est à dire φ^( x ) . Puisqueφ ( (√2 ) ) = ( 2 π)- 1 / deuxe- 1 , votre estimation pour e = ϕ^( 2-√) 2 π--√ .
Si vous voulez devenir totalement fou, vous pouvez même estimer2-√ et 2 π--√ en utilisant les méthodes que vous avez discutées plus tôt.
Méthode nécessitant très peu d'échantillons, mais causant une quantité impie d'erreur numérique
Une réponse complètement stupide, mais très efficace, basée sur un commentaire que j'ai fait:
LaisserX∼ uniforme ( - 1 , 1 ) . DéfinirYn= | ( x¯)n| . Définire^= ( 1 - Yn)- 1 / Yn .
Cela convergera très vite, mais entraînera également une erreur numérique extrême.
Whuber a fait remarquer qu'il utilisait la fonction power, qui appelle généralement la fonction exp. Cela pourrait être évité en discrétisantYn , tel que 1 / Yn is an integer and the power could be replaced with repeated multiplication. It would be required that as n→∞ , the discretizing of Yn would get finer and finer,and the discretization would have to exclude Yn=0 . With this, the estimator theoretically (i.e. the world in which numeric error does not exist) would converge to e , and quite fast!
la source
Voici une autre façon de procéder, même s’il est assez lent. Je ne prétends pas à l'efficacité, mais propose cette alternative dans un esprit de complétude.
Contra Xi'an , je vais supposer, aux fins de cette question, que vous êtes capable de générer et d’utiliser une séquence den variables pseudo-aléatoires uniformes U1, ⋯ , Un~ IID U ( 0 , 1 ) et vous devez ensuite estimer e par une méthode utilisant des opérations arithmétiques de base (en d’autres termes, vous ne pouvez pas utiliser de fonctions logarithmiques ou exponentielles ni de distributions utilisant ces fonctions).† La méthode actuelle est motivée par un résultat simple impliquant des variables aléatoires uniformes:
L'estimatione en utilisant ce résultat: nous ordonnons d'abord les valeurs d'échantillon en ordre décroissant pour obtenir les statistiques d'ordrevous( 1 )⩾ ⋯ ⩾ u( n ) et ensuite nous définissons les sommes partielles:
Maintenant, laissem ≡ min { k | S( k ) ⩾ 1 } puis estimer 1 / e par interpolation des variables uniformes ordonnées. Cela donne un estimateur poure
donné par:
Cette méthode présente un léger biais (dû à l’interpolation linéaire du point de coupure pour1 / e ) mais c’est un estimateur cohérent pour e . La méthode peut être mise en œuvre assez facilement, mais elle nécessite un tri des valeurs, ce qui nécessite davantage de calculs que le calcul déterministe dee . Cette méthode est lente car elle implique le tri des valeurs.
Implémentation dans R: La méthode peut être implémentée en
R
utilisantrunif
pour générer des valeurs uniformes. Le code est comme suit:La mise en œuvre de ce code fait converger la vraie valeur dee , mais il est très lent comparé aux méthodes déterministes.
la source