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 \r
  44 using namespace std;\r
  45 \r
  46 #define Pi 3.1415926535897932\r
  47 \r
  48 double lngamma(double z)\r
  49 {\r
  50         double result,sum;\r
  51         static double c[7]={1.000000000190015, 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5};\r
  52         \r
  53         sum=c[0];\r
  54         for(int i=1;i<=6;i++)sum += c[i]/(z+i);\r
  55         \r
  56         result = (z+0.5)*log(z+5.5)-(z+5.5)+log(2.5066282746310005*sum/z);\r
  57         \r
  58         return result;\r
  59 }\r
  60 \r
  61 \r
  62 \r
  63 int poissonInt(double lambda) \r
  64 {\r
  65         static double mean=-1, waitTime,s,logmean,c;\r
  66         double count, eventTime,ratio,y;\r
  67         \r
  68         static double maxInt=pow(2.0,(double)(8*sizeof(int)-1))-1;\r
  69         \r
  70         if(lambda <0){\r
  71                 cout<<"The mean of Poisson should be >=0, input="<<lambda<<endl;\r
  72                 exit(1);        \r
  73         }\r
  74         else if(lambda < 12){\r
  75                 if(lambda != mean){\r
  76                         mean = lambda;\r
  77                         waitTime = exp(-lambda);        \r
  78                 }\r
  79                 \r
  80                 count = 0;\r
  81                 eventTime = globalRandom.Next();\r
  82                 \r
  83                 while(eventTime > waitTime){\r
  84                         count++;\r
  85                         eventTime *= globalRandom.Next();       \r
  86                 }\r
  87         }\r
  88         else{\r
  89                 if(lambda != mean){\r
  90                         mean = lambda;\r
  91                         s=sqrt(2*lambda);\r
  92                         logmean=log(lambda);\r
  93                         c=lambda*log(lambda)-lngamma(lambda+1);\r
  94                 }\r
  95                 \r
  96                 do{\r
  97                         do{\r
  98                                 y=tan(Pi*globalRandom.Next());\r
  99                                 count = floor(s*y+lambda);\r
 100                         }while(count<0);\r
 101                         \r
 102                         ratio = 0.9*(1+y*y)*exp(count*logmean-lngamma(count+1)-c);\r
 103                         \r
 104                 }while(globalRandom.Next()>ratio);\r
 105                 \r
 106         }\r
 107         \r
 108         if(count>maxInt)count=maxInt;\r
 109         \r
 110         return (int)count;              \r
 111 }\r
 112 \r
 113 \r
 114 \r
 115 \r