GNU Octave  4.0.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
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
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