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