GNU Octave  8.1.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.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 2006-2023 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 /* Original version written by Paul Kienzle distributed as free
27  software in the in the public domain. */
28 
29 /*
30 
31 double randg (a)
32 void fill_randg (a,n,x)
33 
34 Generate a series of standard gamma distributions.
35 
36 See: Marsaglia G and Tsang W (2000), "A simple method for generating
37 gamma variables", ACM Transactions on Mathematical Software 26(3) 363-372
38 
39 Needs the following defines:
40 * NAN: value to return for Not-A-Number
41 * RUNI: uniform generator on (0,1)
42 * RNOR: normal generator
43 * REXP: exponential generator, or -log(RUNI) if one isn't available
44 * INFINITE: function to test whether a value is infinite
45 
46 Test using:
47  mean = a
48  variance = a
49  skewness = 2/sqrt(a)
50  kurtosis = 3 + 6/sqrt(a)
51 
52 Note that randg can be used to generate many distributions:
53 
54 gamma(a,b) for a>0, b>0 (from R)
55  r = b*randg(a)
56 beta(a,b) for a>0, b>0
57  r1 = randg(a,1)
58  r = r1 / (r1 + randg(b,1))
59 Erlang(a,n)
60  r = a*randg(n)
61 chisq(df) for df>0
62  r = 2*randg(df/2)
63 t(df) for 0<df<inf (use randn if df is infinite)
64  r = randn () / sqrt(2*randg(df/2)/df)
65 F(n1,n2) for 0<n1, 0<n2
66  r1 = 2*randg(n1/2)/n1 or 1 if n1 is infinite
67  r2 = 2*randg(n2/2)/n2 or 1 if n2 is infinite
68  r = r1 / r2
69 negative binonial (n, p) for n>0, 0<p<=1
70  r = randp((1-p)/p * randg(n))
71  (from R, citing Devroye(1986), Non-Uniform Random Variate Generation)
72 non-central chisq(df,L), for df>=0 and L>0 (use chisq if L=0)
73  r = randp(L/2)
74  r(r>0) = 2*randg(r(r>0))
75  r(df>0) += 2*randg(df(df>0)/2)
76  (from R, citing formula 29.5b-c in Johnson, Kotz, Balkrishnan(1995))
77 Dirichlet(a1,...,ak) for ai > 0
78  r = (randg(a1),...,randg(ak))
79  r = r / sum(r)
80  (from GSL, citing Law & Kelton(1991), Simulation Modeling and Analysis)
81 */
82 
83 #if defined (HAVE_CONFIG_H)
84 # include "config.h"
85 #endif
86 
87 #include <cmath>
88 
89 #include "lo-ieee.h"
90 #include "randgamma.h"
91 #include "randmtzig.h"
92 
94 
95 template <typename T> void rand_gamma (T a, octave_idx_type n, T *r)
96 {
98  /* If a < 1, start by generating gamma (1+a) */
99  const T d = (a < 1. ? 1.+a : a) - 1./3.;
100  const T c = 1./std::sqrt (9.*d);
101 
102  /* Handle invalid cases */
103  if (a <= 0 || lo_ieee_isinf (a))
104  {
105  for (i=0; i < n; i++)
106  r[i] = numeric_limits<T>::NaN ();
107  return;
108  }
109 
110  for (i=0; i < n; i++)
111  {
112  T x, xsq, v, u;
113  restart:
114  x = rand_normal<T> ();
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 = rand_uniform<T> ();
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 (-rand_exponential<T> () / a);
131  }
132 }
133 
134 template OCTAVE_API void rand_gamma (double, octave_idx_type, double *);
135 template OCTAVE_API void rand_gamma (float, octave_idx_type, float *);
136 
OCTAVE_END_NAMESPACE(octave)
#define NaN
Definition: Faddeeva.cc:261
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
#define lo_ieee_isinf(x)
Definition: lo-ieee.h:108
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
F77_RET_T const F77_DBLE * x
#define OCTAVE_API
Definition: main.in.cc:55
octave_idx_type n
Definition: mx-inlines.cc:753
T * r
Definition: mx-inlines.cc:773
void rand_gamma(T a, octave_idx_type n, T *r)
Definition: randgamma.cc:95