ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
dranpois.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: dranpois.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2009-2012 ADMB Foundation
00006  */
00011 #include "fvar.hpp"
00012 /*
00013 #define IM1 2147483563
00014 #define IM2 2147483399
00015 #define AM (1.0/IM1)
00016 #define IMM1 (IM1-1)
00017 #define IA1 40014
00018 #define IA2 40692
00019 #define IQ1 53668
00020 #define IQ2 52774
00021 #define IR1 12211
00022 #define IR2 3791
00023 #define NTAB 32
00024 #define NDIV (1+IMM1/NTAB)
00025 #define EPS 1.2e-7
00026 #define RNMX (1.0-EPS)
00027 */
00028 
00037 double randpoisson(double xm, const random_number_generator& rng)
00038 {
00039   double gammln(double xx);
00040   static double sq,alxm,g,oldm=(-1.0);
00041 
00042   double em,t,y;
00043 
00044   if (xm < 12.0) {
00045     if (xm != oldm) {
00046       oldm=xm;
00047       g=exp(-xm);
00048     }
00049     em = -1;
00050     t=1.0;
00051     do {
00052       ++em;
00053       //t *= ran1(idum);
00054       t*=  y=((random_number_generator&) rng).better_rand();
00055     } while (t > g);
00056   } else {
00057     if (xm != oldm) {
00058       oldm=xm;
00059       sq=sqrt(2.0*xm);
00060       alxm=log(xm);
00061       g=xm*alxm-gammln(xm+1.0);
00062     }
00063     do {
00064       do {
00065         //y=tan(PI*ran1(idum));
00066         y=tan(PI*((random_number_generator&) rng).better_rand());
00067         em=sq*y+xm;
00068       } while (em < 0.0);
00069       em=floor(em);
00070       t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
00071     } while (((random_number_generator&) rng).better_rand() > t);
00072     //} while (ran1(idum) > t);
00073   }
00074   return em;
00075 }
00076 #undef PI
00077 
00085 void dvector::fill_randpoisson(double lambda,
00086   const random_number_generator& rng)
00087   {
00088     for (int i=indexmin(); i<=indexmax(); i++)
00089     {
00090       elem(i)=randpoisson(lambda,rng);
00091     }
00092   }