|
7019
|
1 /* |
|
|
2 |
|
|
3 Copyright (C) 2006, 2007 John W. Eaton |
|
|
4 |
|
|
5 This file is part of Octave. |
|
|
6 |
|
|
7 Octave is free software; you can redistribute it and/or modify it |
|
|
8 under the terms of the GNU General Public License as published by the |
|
|
9 Free Software Foundation; either version 3 of the License, or (at your |
|
|
10 option) any later version. |
|
|
11 |
|
|
12 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
15 for more details. |
|
|
16 |
|
|
17 You should have received a copy of the GNU General Public License |
|
|
18 along with Octave; see the file COPYING. If not, see |
|
|
19 <http://www.gnu.org/licenses/>. |
|
|
20 |
|
|
21 */ |
|
|
22 |
|
|
23 /* Original version written by Paul Kienzle distributed as free |
|
|
24 software in the in the public domain. */ |
|
5742
|
25 |
|
|
26 /* |
|
|
27 |
|
|
28 double randg(a) |
|
|
29 void fill_randg(a,n,x) |
|
|
30 |
|
|
31 Generate a series of standard gamma distributions. |
|
|
32 |
|
|
33 See: Marsaglia G and Tsang W (2000), "A simple method for generating |
|
|
34 gamma variables", ACM Transactions on Mathematical Software 26(3) 363-372 |
|
|
35 |
|
|
36 Needs the following defines: |
|
|
37 * NAN: value to return for Not-A-Number |
|
|
38 * RUNI: uniform generator on (0,1) |
|
|
39 * RNOR: normal generator |
|
|
40 * REXP: exponential generator, or -log(RUNI) if one isn't available |
|
|
41 * INFINITE: function to test whether a value is infinite |
|
|
42 |
|
|
43 Test using: |
|
|
44 mean = a |
|
|
45 variance = a |
|
|
46 skewness = 2/sqrt(a) |
|
|
47 kurtosis = 3 + 6/sqrt(a) |
|
|
48 |
|
|
49 Note that randg can be used to generate many distributions: |
|
|
50 |
|
|
51 gamma(a,b) for a>0, b>0 (from R) |
|
|
52 r = b*randg(a) |
|
|
53 beta(a,b) for a>0, b>0 |
|
|
54 r1 = randg(a,1) |
|
|
55 r = r1 / (r1 + randg(b,1)) |
|
|
56 Erlang(a,n) |
|
|
57 r = a*randg(n) |
|
|
58 chisq(df) for df>0 |
|
|
59 r = 2*randg(df/2) |
|
|
60 t(df) for 0<df<inf (use randn if df is infinite) |
|
|
61 r = randn() / sqrt(2*randg(df/2)/df) |
|
|
62 F(n1,n2) for 0<n1, 0<n2 |
|
|
63 r1 = 2*randg(n1/2)/n1 or 1 if n1 is infinite |
|
|
64 r2 = 2*randg(n2/2)/n2 or 1 if n2 is infinite |
|
|
65 r = r1 / r2 |
|
|
66 negative binonial (n, p) for n>0, 0<p<=1 |
|
|
67 r = randp((1-p)/p * randg(n)) |
|
|
68 (from R, citing Devroye(1986), Non-Uniform Random Variate Generation) |
|
|
69 non-central chisq(df,L), for df>=0 and L>0 (use chisq if L=0) |
|
|
70 r = randp(L/2) |
|
|
71 r(r>0) = 2*randg(r(r>0)) |
|
|
72 r(df>0) += 2*randg(df(df>0)/2) |
|
|
73 (from R, citing formula 29.5b-c in Johnson, Kotz, Balkrishnan(1995)) |
|
|
74 Dirichlet(a1,...,ak) for ai > 0 |
|
|
75 r = (randg(a1),...,randg(ak)) |
|
|
76 r = r / sum(r) |
|
|
77 (from GSL, citing Law & Kelton(1991), Simulation Modeling and Analysis) |
|
|
78 */ |
|
|
79 |
|
|
80 #if defined (HAVE_CONFIG_H) |
|
|
81 #include <config.h> |
|
|
82 #endif |
|
|
83 |
|
|
84 #include <math.h> |
|
|
85 #include <stdio.h> |
|
|
86 |
|
|
87 #include "lo-ieee.h" |
|
|
88 #include "randmtzig.h" |
|
|
89 #include "randgamma.h" |
|
|
90 |
|
|
91 #undef NAN |
|
|
92 #define NAN octave_NaN |
|
|
93 #define INFINITE lo_ieee_isinf |
|
|
94 #define RUNI oct_randu() |
|
|
95 #define RNOR oct_randn() |
|
|
96 #define REXP oct_rande() |
|
|
97 |
|
|
98 void |
|
|
99 oct_fill_randg (double a, octave_idx_type n, double *r) |
|
|
100 { |
|
|
101 octave_idx_type i; |
|
|
102 /* If a < 1, start by generating gamma(1+a) */ |
|
|
103 const double d = (a < 1. ? 1.+a : a) - 1./3.; |
|
|
104 const double c = 1./sqrt(9.*d); |
|
|
105 |
|
|
106 /* Handle invalid cases */ |
|
|
107 if (a <= 0 || INFINITE(a)) |
|
|
108 { |
|
|
109 for (i=0; i < n; i++) |
|
|
110 r[i] = NAN; |
|
|
111 return; |
|
|
112 } |
|
|
113 |
|
|
114 for (i=0; i < n; i++) |
|
|
115 { |
|
|
116 double x, xsq, v, u; |
|
|
117 restart: |
|
|
118 x = RNOR; |
|
|
119 v = (1+c*x); |
|
|
120 v *= v*v; |
|
|
121 if (v <= 0) |
|
|
122 goto restart; /* rare, so don't bother moving up */ |
|
|
123 u = RUNI; |
|
|
124 xsq = x*x; |
|
|
125 if (u >= 1.-0.0331*xsq*xsq && log(u) >= 0.5*xsq + d*(1-v+log(v))) |
|
|
126 goto restart; |
|
|
127 r[i] = d*v; |
|
|
128 } |
|
|
129 if (a < 1) |
|
|
130 { /* Use gamma(a) = gamma(1+a)*U^(1/a) */ |
|
|
131 /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */ |
|
|
132 for (i = 0; i < n; i++) |
|
|
133 r[i] *= exp(-REXP/a); |
|
|
134 } |
|
|
135 } |
|
|
136 |
|
|
137 double |
|
|
138 oct_randg (double a) |
|
|
139 { |
|
|
140 double ret; |
|
|
141 oct_fill_randg(a,1,&ret); |
|
|
142 return ret; |
|
|
143 } |
|
|
144 |
|
|
145 /* |
|
|
146 ;;; Local Variables: *** |
|
|
147 ;;; mode: C *** |
|
|
148 ;;; End: *** |
|
|
149 */ |
