GNU Octave  4.2.1
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.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2006-2017 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 "lo-ieee.h"
85 #include "lo-math.h"
86 #include "randmtzig.h"
87 #include "randgamma.h"
88 
89 #define INFINITE lo_ieee_isinf
90 #define RUNI oct_randu()
91 #define RNOR oct_randn()
92 #define REXP oct_rande()
93 
94 void
95 oct_fill_randg (double a, octave_idx_type n, double *r)
96 {
98  /* If a < 1, start by generating gamma (1+a) */
99  const double d = (a < 1. ? 1.+a : a) - 1./3.;
100  const double c = 1./sqrt (9.*d);
101 
102  /* Handle invalid cases */
103  if (a <= 0 || INFINITE(a))
104  {
105  for (i=0; i < n; i++)
107  return;
108  }
109 
110  for (i=0; i < n; i++)
111  {
112  double x, xsq, v, u;
113  restart:
114  x = RNOR;
115  v = (1+c*x);
116  v *= v*v;
117  if (v <= 0)
118  goto restart; /* rare, so don't bother moving up */
119  u = RUNI;
120  xsq = x*x;
121  if (u >= 1.-0.0331*xsq*xsq && std::log (u) >= 0.5*xsq + d*(1-v+std::log (v)))
122  goto restart;
123  r[i] = d*v;
124  }
125  if (a < 1)
126  {
127  /* Use gamma(a) = gamma(1+a)*U^(1/a) */
128  /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */
129  for (i = 0; i < n; i++)
130  r[i] *= exp (-REXP/a);
131  }
132 }
133 
134 double
135 oct_randg (double a)
136 {
137  double ret;
138  oct_fill_randg (a,1,&ret);
139  return ret;
140 }
141 
142 #undef RUNI
143 #undef RNOR
144 #undef REXP
145 #define RUNI oct_float_randu()
146 #define RNOR oct_float_randn()
147 #define REXP oct_float_rande()
148 
149 void
151 {
153  /* If a < 1, start by generating gamma(1+a) */
154  const float d = (a < 1. ? 1.+a : a) - 1./3.;
155  const float c = 1./sqrt (9.*d);
156 
157  /* Handle invalid cases */
158  if (a <= 0 || INFINITE(a))
159  {
160  for (i=0; i < n; i++)
162  return;
163  }
164 
165  for (i=0; i < n; i++)
166  {
167  float x, xsq, v, u;
168  frestart:
169  x = RNOR;
170  v = (1+c*x);
171  v *= v*v;
172  if (v <= 0)
173  goto frestart; /* rare, so don't bother moving up */
174  u = RUNI;
175  xsq = x*x;
176  if (u >= 1.-0.0331*xsq*xsq && std::log (u) >= 0.5*xsq + d*(1-v+std::log (v)))
177  goto frestart;
178  r[i] = d*v;
179  }
180  if (a < 1)
181  {
182  /* Use gamma(a) = gamma(1+a)*U^(1/a) */
183  /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */
184  for (i = 0; i < n; i++)
185  r[i] *= exp (-REXP/a);
186  }
187 }
188 
189 float
191 {
192  float ret;
193  oct_fill_float_randg (a,1,&ret);
194  return ret;
195 }
#define RNOR
Definition: randgamma.cc:146
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the base of natural logarithms The constant ex $e satisfies the equation log(e)
void oct_fill_randg(double a, octave_idx_type n, double *r)
Definition: randgamma.cc:95
u
Definition: lu.cc:138
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
#define REXP
Definition: randgamma.cc:147
float oct_float_randg(float a)
Definition: randgamma.cc:190
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
Definition: lu.cc:138
double oct_randg(double a)
Definition: randgamma.cc:135
void oct_fill_float_randg(float a, octave_idx_type n, float *r)
Definition: randgamma.cc:150
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
#define INFINITE
Definition: randgamma.cc:89
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the IEEE symbol NaN(Not a Number).NaN is the result of operations which do not produce a well defined 0 result.Common operations which produce a NaN are arithmetic with infinity ex($\infty-\infty $)
#define RUNI
Definition: randgamma.cc:145
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x