1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2008-2010 The Regents of the University of California
4 //
5 // This file is part of Qbox
6 //
7 // Qbox is distributed under the terms of the GNU General Public License
8 // as published by the Free Software Foundation, either version 2 of
9 // the License, or (at your option) any later version.
10 // See the file COPYING in the root directory of this distribution
11 // or <http://www.gnu.org/licenses/>.
12 //
13 ////////////////////////////////////////////////////////////////////////////////
14 //
15 // sampling.cpp
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18
19 #include <cmath> // sqrt, log
20 #include <cstdlib> // drand48
21 using namespace std;
22
23 ////////////////////////////////////////////////////////////////////////////////
normal_dev(double * d1,double * d2)24 void normal_dev(double* d1, double *d2)
25 {
26 // Box-Muller algorithm: generate a pair of normal random deviates
27 // J. R. Bell, Algorithm 334, Comm. ACM 11, 498, (1968).
28 // modified following R. Knop, Comm. ACM 12, 281 (1969).
29 double x,y,s;
30 do
31 {
32 x = 2.0 * drand48() - 1.0;
33 y = 2.0 * drand48() - 1.0;
34 s = x*x + y*y;
35 }
36 while ( s > 1.0 );
37 double l = sqrt(-2.0 * log(s)/s);
38 *d1 = x * l;
39 *d2 = y * l;
40 }
41
42 ////////////////////////////////////////////////////////////////////////////////
43 // generate deviates for the gamma function
44 // use the sum of nsum gaussian deviates
45 // use: tgamma_sum nsum nval seed
46 ////////////////////////////////////////////////////////////////////////////////
gamma_dev(int a)47 double gamma_dev(int a)
48 {
49 double am=a-1.0, e, s=sqrt(2.0*(a-1.0)+1.0), v1, v2, x, y;
50 do
51 {
52 do
53 {
54 do
55 {
56 v1 = 2.0*drand48()-1.0;
57 v2 = 2.0*drand48()-1.0;
58 } while ( v1*v1 + v2*v2 > 1.0 || v1 == 0.0 );
59 y = v2/v1;
60 x = s*y+am;
61 } while ( x <= 0.0 );
62 e = (1.0+y*y)*exp(am*log(x/am)-s*y);
63 } while ( drand48() > e );
64 return x;
65 }
66