1 //////////////////////////////////////////////////////////////////////
\r
3 //////////////////////////////////////////////////////////////////////////////
\r
4 // COPYRIGHT NOTICE FOR GENOME CODE
\r
6 // Copyright (C) 2006 - 2009, Liming Liang and Goncalo Abecasis,
\r
7 // All rights reserved.
\r
9 // Redistribution and use in source and binary forms, with or without
\r
10 // modification, are permitted provided that the following conditions
\r
13 // 1. Redistributions of source code must retain the above copyright
\r
14 // notice, this list of conditions and the following disclaimer.
\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
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
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
36 ///////////////////////////////////////////////////////////////////////////////
\r
45 using namespace std;
\r
47 #define Pi 3.1415926535897932
\r
49 double lngamma(double z)
\r
52 static double c[7]={1.000000000190015, 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5};
\r
55 for(int i=1;i<=6;i++)sum += c[i]/(z+i);
\r
57 result = (z+0.5)*log(z+5.5)-(z+5.5)+log(2.5066282746310005*sum/z);
\r
64 int poissonInt(double lambda)
\r
66 static double mean=-1, waitTime,s,logmean,c;
\r
67 double count, eventTime,ratio,y;
\r
69 static double maxInt=pow(2.0,(double)(8*sizeof(int)-1))-1;
\r
72 cout<<"The mean of Poisson should be >=0, input="<<lambda<<endl;
\r
75 else if(lambda < 12){
\r
78 waitTime = exp(-lambda);
\r
82 eventTime = globalRandom.Next();
\r
84 while(eventTime > waitTime){
\r
86 eventTime *= globalRandom.Next();
\r
93 logmean=log(lambda);
\r
94 c=lambda*log(lambda)-lngamma(lambda+1);
\r
99 y=tan(Pi*globalRandom.Next());
\r
100 count = floor(s*y+lambda);
\r
103 ratio = 0.9*(1+y*y)*exp(count*logmean-lngamma(count+1)-c);
\r
105 }while(globalRandom.Next()>ratio);
\r
109 if(count>maxInt)count=maxInt;
\r
111 return (int)count;
\r