## Гнеде́нко and Forsythe [and e]

**I**n the wake of my earlier post on the Monte Carlo estimation of *e* and *e⁻¹*, after a discussion with my colleague Murray Pollock (Warwick) Gnedenko’s solution, he pointed out another (early) Monte Carlo approximation called Forsythe’s method. That is detailed quite neatly in Luc Devroye’s bible, Non-uniform random variate generation (a free bible!). The idea is to run a sequence of uniform generations until the sequence goes up, i.e., until the latest uniform is larger than the previous one. The expectation of the corresponding stopping rule, N, which is the number of generations the uniform sequence went down continuously is then *e*, while the probability that N is odd is *e⁻¹, *most unexpectedly! Forsythe’s method actually aims at a much broader goal, namely simulating from any density of the form g(x) exp{-F(x)}, F taking values in (0,1). This method generalises von Neuman’s exponential generator (see Devroye, p.126) which only requires uniform generations.

This is absolutely identical to Gnedenko’s approach in that both events have the same 1/n! probability to occur [as pointed out by Gérard Letac in a comment on the previous entry]. *(I certainly cannot say whether or not one of the authors was aware of the other’s result: Forsythe generalised von Neumann‘s method around 1972, while Gnedenko published Theory of Probability at least in 1969, but this may be the date of the English translation, I have not been able to find the reference on the Russian wikipedia page…)* Running a small R experiment to compare both distributions of N, the above barplot shows that they look quite similar:

n=1e6 use=runif(n) # Gnedenko's in action: gest=NULL i=1 while (i<(n-100)){ sumuse=cumsum(use[i:(i+10)]) if (sumuse[11]<1]) sumuse=cumsum(use[i:(i+100)]) j=min((1:length(sumuse))[sumuse>1]) gest=c(gest,j) i=i+j} #Forsythe's method: fest=NULL i=1 while (i<(n-100)){ sumuse=c(-1,diff(use[i:(i+10)])) if (max(sumuse)<0]) sumuse=c(-1,diff(use[i:(i+100)])) j=min((1:length(sumuse))[sumuse>0]) fest=c(fest,j) i=i+j}

And the execution times of both approaches [with this rudimentary R code!] are quite close.

February 23, 2016 at 3:54 pm

I haven’t run a speed test, but this code seems a lot simpler. I ran some cheks, and never got a decreasing run length greater than 10.

n=25

fest=NULL

for (k in 1:10000){

use=runif(n)

fest[k] 0)[1] + 1

}

mean(fest)

length(fest)/sum(fest%%2)

February 23, 2016 at 3:55 pm

sorry, guess a “\0)[1] + 1

}

February 23, 2016 at 3:56 pm

still can’t get my code to paste in correctly . try again:

n=25

fest=NULL

for (k in 1:10000){

use=runif(n)

fest[k] = which(diff(use)>0)[1] + 1

}

February 23, 2016 at 4:49 pm

Thanks Carl for the faster R code. It took me a wee while to spot the [1] after the which and to agree this R code produces the right output. It looks rather similar to the improved code I posted yesterday, I believe.

February 19, 2016 at 12:09 am

I encountered Gnedenko’s book on probability over 30 years ago. A classic with some wonderful stuff that’s not in your standard treatments–especially today.

February 16, 2016 at 7:33 pm

In my copy of the translation of Gnedenko’s 4th edition this is problem 22 on p 235. Our library catalogue has Russian editions back 1950, but whether earlier editions contain the problem, I don’t know.