approximatif en

35

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 e [nombre d'Euler] en utilisant l'intégration de Monte Carlo.

Avez-vous des indications sur la façon dont cela peut être fait?

statisticnewbie12345
la source
7
Il y a beaucoup, beaucoup de façons de faire cela. Cela pourrait devenir évident en contemplant ce que le Rcommandement 2 + mean(exp(-lgamma(ceiling(1/runif(1e5))-1)))fait. (Si vous utilisez la fonction de consignation Gamma vous dérange, remplacez-la par 2 + 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 estimer e avec une précision donnée?
whuber
4
Quelle question délicieuse! J'ai hâte de lire les réponses des autres. Une façon d'attirer l'attention sur cette question - peut-être une demi-douzaine de réponses - serait de réviser la question et de demander des réponses efficaces , comme suggéré par Whuber. C'est comme l'herbe à chat pour les utilisateurs de CV.
Sycorax dit de réintégrer Monica le
1
@EngrStudent Je ne suis pas sûr que l'analogue géométrique existe pour e . Ce n'est tout simplement pas une quantité géométrique naturelle (jeu de mots) comme π .
Aksakal
6
@Aksakal e est une quantité exceptionnellement géométrique. Au niveau le plus élémentaire, il apparaît naturellement dans les expressions pour les domaines liés aux hyperboles. À un niveau légèrement plus avancé, il est intimement lié à toutes les fonctions périodiques, y compris les fonctions trigonométriques, dont le contenu géométrique est évident. Le véritable défi ici est qu’il est si facile de simuler des valeurs liées à e !
whuber
2
@StatsStudent: e n'est pas intéressant en soi. Cependant, si cela conduit à des estimateurs sans biais de quantités comme
exp{0xf(y)dG(y)}
ce qui peut se révéler très utile pour les algorithmes MCMC.
Xi'an

Réponses:

34

La méthode simple et élégante pour estimer e 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 que E[ξ]=e , où ξ est une variable aléatoire qui est définie comme suit. C'est le nombre minimum de n tel que i=1nri>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:

 0.0180
 0.4596
 0.7920

Ensuite, votre premier essai rendu 3. Continuez à faire ces essais, et vous remarquerez qu'en moyenne, vous obtenez e .

Le code MATLAB, le résultat de la simulation et l'histogramme suivent.

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)

Le résultat et l'histogramme:

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

entrez la description de l'image ici

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:

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

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:

entrez la description de l'image ici

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 .

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

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:

entrez la description de l'image ici

É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.

Aksakal
la source
Wow, c'est cool! Pourriez-vous ajouter un paragraphe ou deux expliquant pourquoi cela fonctionne?
Sycorax dit Rétablir Monica
7
(+1) brillant! La réponse mérite la note la plus élevée, car elle ne repose que sur des simulations uniformes. Et n'utilise aucune approximation mais celle due à Monte Carlo. Le fait qu'il soit connecté à Gnedenko est un avantage supplémentaire.
Xi'an
2
Cool! Voici le code Mathematica pour idem, comme un one-liner:
Mean[Table[ Length[NestWhileList[(Random[]+#) &, Random[], #<1&]], {10^6}]]
wolfies
4
@wolfies La traduction directe suivante de la Rsolution 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]]
whuber
1
J'ai posté le "pourquoi est la moyenne ?" question comme une question à part entière ; Je soupçonne que ma solution de croquis (qui est ce qui m'est immédiatement venu à l'esprit comme une visualisation "évidente" du problème) n'est pas nécessairement la façon dont les étudiants russes étaient censés le faire! Des solutions alternatives seraient donc les bienvenues. e
Silverfish
19

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 de e .

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 .loge

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 - 1nni1n1np=(11n)n.

Parce que je suppose que nous savons que

exp(1)=limn(11n)n

donc on peut aussi écrire

exp(1)p^=i=1mIiBjm

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.pmBji

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, pnp^ 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 - 1e est précis à la troisième décimale.nppnnen=103p1e

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 . mp^pep^entrez la description de l'image ici

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 $.

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 dit Réintégrer Monica
la source
4
+1 Cela fait beaucoup de sens. Avez-vous une chance de partager votre code si vous l'avez écrit?
Antoni Parellada
2
Même si cela peut être arbitrairement exact, il n’est finalement pas satisfaisant car il ne simule qu’une approximation de plutôt que de e lui-même. ee
whuber
1
Sûr. Vous voudriez juste terminer avec un appel répliqué à l'intérieur d'un autre, qui est essentiellement le même que nous avons maintenant.
Sycorax dit: Réintégrer Monica le
1
@whuber Je ne vois pas vraiment la distinction entre une approximation arbitrairement exacte d'une approximation arbitrairement exacte de et une approximation arbitrairement exacte de e elle-même. ee
Jwg
1
@jwg En plus d'être important sur le plan conceptuel, il est également important sur le plan pratique car la mise en œuvre d'une approximation approximative nécessite de garder une trace de la précision de chacune des deux approximations. Mais je conviens que lorsque les deux approximations sont acceptables, l’approche globale est satisfaisante.
whuber
14

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).

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

Remarque 1: Comme indiqué dans les commentaires, il s'agit d'un argument plutôt compliqué, car simuler à partir d'une distribution de Poisson ou, de manière équivalente, d'une distribution exponentielle peut être difficile à imaginer sans impliquer une fonction log ou une fonction exp ... Mais ensuite, Hu Hu sauvetage de cette réponse avec une solution des plus élégantes basée sur des uniformes ordonnés. Ce qui est une approximation cependant, puisque la distribution d'un espacement uniforme est un Beta B ( 1 , n )U(i:n)U(i1:n)B(1,n), ce qui implique que qui converge verse-1lorsquendevient infini. De plus, le générateur exponentiel de von Neumann (1951)utiliseuniquementdes générations uniformes.

P(n{U(je:n)-U(je-1:n)}1)=(1-1n)n
e-1n

Solution 2:

e

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

Solution 3:

u1,u2,... until un+1>un. The expectation of the corresponding stopping rule, N, which is the number of time the uniform sequence went down is then e while the probability that N is odd is e1! (Forsythe's method actually aims at simulating from any density of the form expG(x), hence is more general than approximating e and e1.)

This is quite parallel to Gnedenko's approach used in Aksakal's answer, so I wonder if one can be derived from the other. At the very least, both have the same distribution with probability mass 1/n! for value n.

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:

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
la source
12
As long as one knows how to do Poisson simulation without knowing e.
Glen_b -Reinstate Monica
5
If I call R rpoiss() generator, I can pretend I do not know e. More seriously, you can generate exponential variates E(1) [using a log function rather than e] until the sum exceeds 1 and the resulting number minus one is a Poisson P(1).
Xi'an
5
Computing log is tantamount to computing exp, since they are inverses. You can avoid computing any such function in various ways. Here is one solution based on your first answer: n <- 1e5; 1/mean(n*diff(sort(runif(n+1))) > 1) It uses only elementary arithmetic.
whuber
3
Je crois que la méthode de Forsythe est la même que celle de Gnedenko. Choisir un uniformeXn tel que ΣnXje est inférieur à 1 revient à choisir Xn plus petit que 1-Σn-1Xjeet si nous réussissons, 1-ΣnXje est conditionnellement uniformément réparti entre 1-Σn-1Xje et 0.
jwg
3
I wasn't aware of Forsythe's approach. However, it's linked to something else very interesting. If instead of stopping at n+1 you keep sampling, then the expectation of the distance from n to the next bottom is exactly 3.
Aksakal
7

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:

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

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

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:

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

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:

  • La méthode d'Aksakal: dataA
  • Méthode 1 de Xian utilisant le code whuber: dataB

... both in Mathematica. The following diagram contrasts a nonparametric kernel density estimate of the resulting dataA and dataB data sets.

enter image description here

So, while whuber's code (red curve) is about twice as fast, the method does not appear to be as 'reliable'.

wolfies
la source
Une ligne verticale à l'emplacement de la valeur réelle améliorerait considérablement cette image.
Sycorax dit: Réintégrer Monica le
1
C'est une observation très intéressante, merci. Étant donné que la demi-largeur sera mise à l'échelle de manière quadratique avec la taille de la simulation et que la demi-largeur de la méthode de Xi'an est environ deux fois celle de la méthode d'Aksakal, le fait de parcourir quatre fois plus d'itérations les rend tout aussi précises. La question de savoir combien d'effort est nécessaire à chaque itération demeure: si une itération de la méthode de Xi'an prend moins d'un quart de l'effort, cette méthode serait toujours plus efficace.
whuber
1
Je crois que la situation devient claire lorsque vous comparez le nombre de réalisations de variables aléatoires requises dans les deux méthodes plutôt que la valeur nominale de n.
whuber
1
@whuber a écrit: 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 dedix6 à 6 x dix6 (soit 6 fois le nombre de points) produit une courbe similaire à celle d’Aksaksal.
loups
1
Bien fait avec le code - il sera difficile d'améliorer beaucoup à ce sujet.
whuber
2

Mé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)=eXou 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¯12n~˙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é de N(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.eXrè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'un N(0,1), c'est à dire φ^(X). Puisqueφ((2))=(2π)-1/2e-1, votre estimation pour e=φ^(2)2π.

Si vous voulez devenir totalement fou, vous pouvez même estimer 2 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:

Laisser X~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!

Cliff AB
la source
2
L’approche CLT n’est pas satisfaisante, car vous savez que ces valeurs ne sont pas distribuées normalement. Mais il existe de nombreuses façons de générer des variables normales sans avoir besoin deeou logarithmes: la méthode de Box-Muller en est une. Celui-ci, cependant, nécessite des fonctions trigonométriques et (à un niveau fondamental) identiques à des exponentielles.
whuber
1
@whuber: Je n'ai pas utilisé Box-Muller en raison de la transformation du journal requise trop directement en exponentielle dans mon livre. Je me suis permis réflexivement cos et le péché, mais ce fut seulement parce que j'avais oublié l' analyse complexe pour un moment, donc bon point.
Cliff AB
1
However, I would take argument with the idea that the generated normal approximation is the weak point of this idea; the density estimation is even weaker! You can think of this idea of having two parameters: n1, the number uniforms used in your "approximated normal" and n2 the number of approximated normals used estimate the density at ϕ(2). As both n1 and n2 approach , the estimator will approach e. In fact, I'm very confident the convergence rate would be much more limited by n2 than n1; non-parametric density has a slow convergence rate!
Cliff AB
2

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:

E(je(Uje1/e)Uje)=1/e1vousvous=1.

L'estimation een utilisant ce résultat: nous ordonnons d'abord les valeurs d'échantillon en ordre décroissant pour obtenir les statistiques d'ordrevous(1)vous(n) et ensuite nous définissons les sommes partielles:

Sn(k)1nΣje=1k1vous(je)pour tous k=1,..,n.

Maintenant, laisse mmin{k|S(k)1} puis estimer 1/epar interpolation des variables uniformes ordonnées. Cela donne un estimateur poure donné par:

e^2vous(m)+vous(m+1).

Cette méthode présente un léger biais (dû à l’interpolation linéaire du point de coupure pour 1/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 Rutilisant runifpour générer des valeurs uniformes. Le code est comme suit:

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 mise en œuvre de ce code fait converger la vraie valeur de e, mais il est très lent comparé aux méthodes déterministes.

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

J’estime que nous voulons éviter toute méthode qui utilise une transformation impliquant un exponentiel ou un logarithme. Si nous pouvons utiliser des densités qui utilisent des exponentielles dans leur définition, il est possible de déduiree à partir de ceux-ci algébriquement en utilisant un appel de densité.

Rétablir Monica
la source