|
5742
|
1 /* This code is in the public domain */ |
|
|
2 |
|
|
3 /* |
|
|
4 |
|
|
5 double randg(a) |
|
|
6 void fill_randg(a,n,x) |
|
|
7 |
|
|
8 Generate a series of standard gamma distributions. |
|
|
9 |
|
|
10 See: Marsaglia G and Tsang W (2000), "A simple method for generating |
|
|
11 gamma variables", ACM Transactions on Mathematical Software 26(3) 363-372 |
|
|
12 |
|
|
13 Needs the following defines: |
|
|
14 * NAN: value to return for Not-A-Number |
|
|
15 * RUNI: uniform generator on (0,1) |
|
|
16 * RNOR: normal generator |
|
|
17 * REXP: exponential generator, or -log(RUNI) if one isn't available |
|
|
18 * INFINITE: function to test whether a value is infinite |
|
|
19 |
|
|
20 Test using: |
|
|
21 mean = a |
|
|
22 variance = a |
|
|
23 skewness = 2/sqrt(a) |
|
|
24 kurtosis = 3 + 6/sqrt(a) |
|
|
25 |
|
|
26 Note that randg can be used to generate many distributions: |
|
|
27 |
|
|
28 gamma(a,b) for a>0, b>0 (from R) |
|
|
29 r = b*randg(a) |
|
|
30 beta(a,b) for a>0, b>0 |
|
|
31 r1 = randg(a,1) |
|
|
32 r = r1 / (r1 + randg(b,1)) |
|
|
33 Erlang(a,n) |
|
|
34 r = a*randg(n) |
|
|
35 chisq(df) for df>0 |
|
|
36 r = 2*randg(df/2) |
|
|
37 t(df) for 0<df<inf (use randn if df is infinite) |
|
|
38 r = randn() / sqrt(2*randg(df/2)/df) |
|
|
39 F(n1,n2) for 0<n1, 0<n2 |
|
|
40 r1 = 2*randg(n1/2)/n1 or 1 if n1 is infinite |
|
|
41 r2 = 2*randg(n2/2)/n2 or 1 if n2 is infinite |
|
|
42 r = r1 / r2 |
|
|
43 negative binonial (n, p) for n>0, 0<p<=1 |
|
|
44 r = randp((1-p)/p * randg(n)) |
|
|
45 (from R, citing Devroye(1986), Non-Uniform Random Variate Generation) |
|
|
46 non-central chisq(df,L), for df>=0 and L>0 (use chisq if L=0) |
|
|
47 r = randp(L/2) |
|
|
48 r(r>0) = 2*randg(r(r>0)) |
|
|
49 r(df>0) += 2*randg(df(df>0)/2) |
|
|
50 (from R, citing formula 29.5b-c in Johnson, Kotz, Balkrishnan(1995)) |
|
|
51 Dirichlet(a1,...,ak) for ai > 0 |
|
|
52 r = (randg(a1),...,randg(ak)) |
|
|
53 r = r / sum(r) |
|
|
54 (from GSL, citing Law & Kelton(1991), Simulation Modeling and Analysis) |
|
|
55 */ |
|
|
56 |
|
|
57 #if defined (HAVE_CONFIG_H) |
|
|
58 #include <config.h> |
|
|
59 #endif |
|
|
60 |
|
|
61 #include <math.h> |
|
|
62 #include <stdio.h> |
|
|
63 |
|
|
64 #include "lo-ieee.h" |
|
|
65 #include "randmtzig.h" |
|
|
66 #include "randgamma.h" |
|
|
67 |
|
|
68 #undef NAN |
|
|
69 #define NAN octave_NaN |
|
|
70 #define INFINITE lo_ieee_isinf |
|
|
71 #define RUNI oct_randu() |
|
|
72 #define RNOR oct_randn() |
|
|
73 #define REXP oct_rande() |
|
|
74 |
|
|
75 void |
|
|
76 oct_fill_randg (double a, octave_idx_type n, double *r) |
|
|
77 { |
|
|
78 octave_idx_type i; |
|
|
79 /* If a < 1, start by generating gamma(1+a) */ |
|
|
80 const double d = (a < 1. ? 1.+a : a) - 1./3.; |
|
|
81 const double c = 1./sqrt(9.*d); |
|
|
82 |
|
|
83 /* Handle invalid cases */ |
|
|
84 if (a <= 0 || INFINITE(a)) |
|
|
85 { |
|
|
86 for (i=0; i < n; i++) |
|
|
87 r[i] = NAN; |
|
|
88 return; |
|
|
89 } |
|
|
90 |
|
|
91 for (i=0; i < n; i++) |
|
|
92 { |
|
|
93 double x, xsq, v, u; |
|
|
94 restart: |
|
|
95 x = RNOR; |
|
|
96 v = (1+c*x); |
|
|
97 v *= v*v; |
|
|
98 if (v <= 0) |
|
|
99 goto restart; /* rare, so don't bother moving up */ |
|
|
100 u = RUNI; |
|
|
101 xsq = x*x; |
|
|
102 if (u >= 1.-0.0331*xsq*xsq && log(u) >= 0.5*xsq + d*(1-v+log(v))) |
|
|
103 goto restart; |
|
|
104 r[i] = d*v; |
|
|
105 } |
|
|
106 if (a < 1) |
|
|
107 { /* Use gamma(a) = gamma(1+a)*U^(1/a) */ |
|
|
108 /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */ |
|
|
109 for (i = 0; i < n; i++) |
|
|
110 r[i] *= exp(-REXP/a); |
|
|
111 } |
|
|
112 } |
|
|
113 |
|
|
114 double |
|
|
115 oct_randg (double a) |
|
|
116 { |
|
|
117 double ret; |
|
|
118 oct_fill_randg(a,1,&ret); |
|
|
119 return ret; |
|
|
120 } |
|
|
121 |
|
|
122 /* |
|
|
123 ;;; Local Variables: *** |
|
|
124 ;;; mode: C *** |
|
|
125 ;;; End: *** |
|
|
126 */ |
