GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
randgamma.c
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2006-2015 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. */
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 <stdio.h>
85 
86 #include "lo-ieee.h"
87 #include "lo-math.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  {
131  /* Use gamma(a) = gamma(1+a)*U^(1/a) */
132  /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */
133  for (i = 0; i < n; i++)
134  r[i] *= exp (-REXP/a);
135  }
136 }
137 
138 double
139 oct_randg (double a)
140 {
141  double ret;
142  oct_fill_randg (a,1,&ret);
143  return ret;
144 }
145 
146 #undef NAN
147 #undef RUNI
148 #undef RNOR
149 #undef REXP
150 #define NAN octave_Float_NaN
151 #define RUNI oct_float_randu()
152 #define RNOR oct_float_randn()
153 #define REXP oct_float_rande()
154 
155 void
156 oct_fill_float_randg (float a, octave_idx_type n, float *r)
157 {
158  octave_idx_type i;
159  /* If a < 1, start by generating gamma(1+a) */
160  const float d = (a < 1. ? 1.+a : a) - 1./3.;
161  const float c = 1./sqrt (9.*d);
162 
163  /* Handle invalid cases */
164  if (a <= 0 || INFINITE(a))
165  {
166  for (i=0; i < n; i++)
167  r[i] = NAN;
168  return;
169  }
170 
171  for (i=0; i < n; i++)
172  {
173  float x, xsq, v, u;
174  frestart:
175  x = RNOR;
176  v = (1+c*x);
177  v *= v*v;
178  if (v <= 0)
179  goto frestart; /* rare, so don't bother moving up */
180  u = RUNI;
181  xsq = x*x;
182  if (u >= 1.-0.0331*xsq*xsq && log (u) >= 0.5*xsq + d*(1-v+log (v)))
183  goto frestart;
184  r[i] = d*v;
185  }
186  if (a < 1)
187  {
188  /* Use gamma(a) = gamma(1+a)*U^(1/a) */
189  /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */
190  for (i = 0; i < n; i++)
191  r[i] *= exp (-REXP/a);
192  }
193 }
194 
195 float
197 {
198  float ret;
199  oct_fill_float_randg (a,1,&ret);
200  return ret;
201 }
#define RUNI
Definition: randgamma.c:151
F77_RET_T const double const double double * d
float oct_float_randg(float a)
Definition: randgamma.c:196
#define REXP
Definition: randgamma.c:153
#define INFINITE
Definition: randgamma.c:93
void oct_fill_float_randg(float a, octave_idx_type n, float *r)
Definition: randgamma.c:156
double oct_randg(double a)
Definition: randgamma.c:139
#define NAN
Definition: randgamma.c:150
#define RNOR
Definition: randgamma.c:152
void oct_fill_randg(double a, octave_idx_type n, double *r)
Definition: randgamma.c:99
F77_RET_T const double * x