The hint’s in the title. I wrote a simple function in Fortran for simulating (or sampling) Poisson random variables. (More precisely, I should say that the function generates Poisson variates.) I used the simple direct method. This method is based on the exponential inter-arrival times of the Poisson (stochastic) process.
My code should not be used for large Poisson parameter values (larger than, say, 20 or 30), as the code may be too slow. Other methods exist for larger parameter values, which I’ve discussed previously.
I just use the standard Fortran function random_number for generating (pseudo-)random numbers. I am not an expert in Fortran, but my Poisson function seems to work fine. I wrote and ran a simple test that estimates the first and second moments, which should match for Poisson variables.
My Fortran code is very similar to the code that I wrote in C and C#, which is located here. You should be able to run it on this website or similar ones that can compile Fortran (95) code.
For various Poisson simulation methods, see the stochastic simulation books by Devroye (Section X.3) or Fishman (Section 8.16). There’s a free online version of Devroye’s book here. The book by Gentle (Section 5.2.8) also briefly covers Poisson variables.
In this post on generating Poisson variates, John D. Cook briefly discusses the direct method (for small Poisson parameter values), as well as a rejection method from a 1979 paper by Atkinson.
I wrote the Poisson code using Fortran 95. There are various books and websites on Fortran. The website tutorialspoint.com gives a good introduction to Fortran. You can also edit, compile and run your Fortran code there with its online Fortran editor. I found this short summary a good start. For alternative Fortran code of a Poisson generator, consult the classic book Numerical Recipes, though I believe the book versions only exist for Fortran 77 and Fortran 90.
On this page I’ve only included the code of the functions for generating uniform and Poisson variates. The rest of the code, including the test, is located here.
!Poisson function -- returns a single Poisson random variable function funPoissonSingle(lambda) result(randPoisson) real, intent(in) :: lambda !input real :: exp_lambda !constant for terminating loop real :: randUni !uniform variable real :: prodUni !product of uniform variables integer :: randPoisson !Poisson variable exp_lambda= exp(-lambda) !initialize variables randPoisson = -1; prodUni = 1; do while (prodUni > exp_lambda) randUni = funUniformSingle() !generate uniform variable prodUni = prodUni * randUni; !update product randPoisson = randPoisson + 1 !increase Poisson variable end do end function !Uniform function -- returns a standard uniform random variable function funUniformSingle() result(randUni) real randUni; call random_seed call random_number(randUni) end function
5 thoughts on “Simulating Poisson random variables in Fortran”
Thanks for the code. Semicolons at the end of a line are redundant. Fortran is not C.
Thanks. Old habits die hard. I do the same with Python, though I read I shouldn’t.
Thanks for the article and code. I was wondering if you have extended the code for large lambda in order of 10^30?
No. I haven’t. You wouldn’t use the direct approach implemented here. You would use one of the approximation techniques that I’ve mentioned in this post:
Check if the older papers have any Fortan code. Happy coding.
Thanks for the reply. I’ve managed to write it based on Ahrens and Dieter’s paper. It took less time than I expected and the algorithm is quite fast!