stochastic.cpp

   1 ////////////////////////////////////////////////////////////////////// \r
   2 // stochastic.cpp \r
   3 //////////////////////////////////////////////////////////////////////////////\r
   4 //              COPYRIGHT NOTICE FOR GENOME CODE\r
   5 //\r
   6 // Copyright (C) 2006 - 2009, Liming Liang and Goncalo Abecasis,\r
   7 // All rights reserved.\r
   8 //\r
   9 //   Redistribution and use in source and binary forms, with or without\r
  10 //   modification, are permitted provided that the following conditions\r
  11 //   are met:\r
  12 //\r
  13 //     1. Redistributions of source code must retain the above copyright\r
  14 //        notice, this list of conditions and the following disclaimer.\r
  15 //\r
  16 //     2. Redistributions in binary form must reproduce the above copyright\r
  17 //        notice, this list of conditions and the following disclaimer in the\r
  18 //        documentation and/or other materials provided with the distribution.\r
  19 //\r
  20 //     3. The names of its contributors may not be used to endorse or promote\r
  21 //        products derived from this software without specific prior written\r
  22 //        permission.\r
  23 //\r
  24 //   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS\r
  25 //   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT\r
  26 //   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR\r
  27 //   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR\r
  28 //   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,\r
  29 //   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,\r
  30 //   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR\r
  31 //   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF\r
  32 //   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING\r
  33 //   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS\r
  34 //   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.\r
  35 //\r
  36 ///////////////////////////////////////////////////////////////////////////////\r
  37 \r
  38 \r
  39 \r
  40 #include "Random.h"\r
  41 #include <math.h>\r
  42 #include <iostream>\r
  43 #include <stdlib.h>\r
  44 \r
  45 using namespace std;\r
  46 \r
  47 #define Pi 3.1415926535897932\r
  48 \r
  49 double lngamma(double z)\r
  50 {\r
  51         double result,sum;\r
  52         static double c[7]={1.000000000190015, 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5};\r
  53         \r
  54         sum=c[0];\r
  55         for(int i=1;i<=6;i++)sum += c[i]/(z+i);\r
  56         \r
  57         result = (z+0.5)*log(z+5.5)-(z+5.5)+log(2.5066282746310005*sum/z);\r
  58         \r
  59         return result;\r
  60 }\r
  61 \r
  62 \r
  63 \r
  64 int poissonInt(double lambda) \r
  65 {\r
  66         static double mean=-1, waitTime,s,logmean,c;\r
  67         double count, eventTime,ratio,y;\r
  68         \r
  69         static double maxInt=pow(2.0,(double)(8*sizeof(int)-1))-1;\r
  70         \r
  71         if(lambda <0){\r
  72                 cout<<"The mean of Poisson should be >=0, input="<<lambda<<endl;\r
  73                 exit(1);        \r
  74         }\r
  75         else if(lambda < 12){\r
  76                 if(lambda != mean){\r
  77                         mean = lambda;\r
  78                         waitTime = exp(-lambda);        \r
  79                 }\r
  80                 \r
  81                 count = 0;\r
  82                 eventTime = globalRandom.Next();\r
  83                 \r
  84                 while(eventTime > waitTime){\r
  85                         count++;\r
  86                         eventTime *= globalRandom.Next();       \r
  87                 }\r
  88         }\r
  89         else{\r
  90                 if(lambda != mean){\r
  91                         mean = lambda;\r
  92                         s=sqrt(2*lambda);\r
  93                         logmean=log(lambda);\r
  94                         c=lambda*log(lambda)-lngamma(lambda+1);\r
  95                 }\r
  96                 \r
  97                 do{\r
  98                         do{\r
  99                                 y=tan(Pi*globalRandom.Next());\r
 100                                 count = floor(s*y+lambda);\r
 101                         }while(count<0);\r
 102                         \r
 103                         ratio = 0.9*(1+y*y)*exp(count*logmean-lngamma(count+1)-c);\r
 104                         \r
 105                 }while(globalRandom.Next()>ratio);\r
 106                 \r
 107         }\r
 108         \r
 109         if(count>maxInt)count=maxInt;\r
 110         \r
 111         return (int)count;              \r
 112 }\r
 113 \r
 114 \r
 115 \r
 116 \r