Mercurial > octave
annotate liboctave/numeric/randgamma.cc @ 29359:7854d5752dd2
maint: merge stable to default.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 10 Feb 2021 10:10:40 -0500 |
parents | 5c14f81e0937 0a5b15007766 |
children | 796f54d4ddbf |
rev | line source |
---|---|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
1 //////////////////////////////////////////////////////////////////////// |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
2 // |
29358
0a5b15007766
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
27923
diff
changeset
|
3 // Copyright (C) 2006-2021 The Octave Project Developers |
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
4 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
5 // See the file COPYRIGHT.md in the top-level directory of this |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
6 // distribution or <https://octave.org/copyright/>. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
7 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
8 // This file is part of Octave. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
9 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
10 // Octave is free software: you can redistribute it and/or modify it |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
11 // under the terms of the GNU General Public License as published by |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
12 // the Free Software Foundation, either version 3 of the License, or |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
13 // (at your option) any later version. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
14 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
15 // Octave is distributed in the hope that it will be useful, but |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
16 // WITHOUT ANY WARRANTY; without even the implied warranty of |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
18 // GNU General Public License for more details. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
19 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
20 // You should have received a copy of the GNU General Public License |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
21 // along with Octave; see the file COPYING. If not, see |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
22 // <https://www.gnu.org/licenses/>. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
23 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
24 //////////////////////////////////////////////////////////////////////// |
7019 | 25 |
26 /* Original version written by Paul Kienzle distributed as free | |
27 software in the in the public domain. */ | |
5742 | 28 |
29 /* | |
30 | |
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
31 double randg (a) |
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
32 void fill_randg (a,n,x) |
5742 | 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 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
39 Needs the following defines: |
5742 | 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) | |
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
64 r = randn () / sqrt(2*randg(df/2)/df) |
5742 | 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) | |
21301
40de9f8f23a6
Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents:
21202
diff
changeset
|
84 # include "config.h" |
5742 | 85 #endif |
86 | |
23662
bd77ab816e43
eliminate obsolete file lo-math.h
John W. Eaton <jwe@octave.org>
parents:
23649
diff
changeset
|
87 #include <cmath> |
bd77ab816e43
eliminate obsolete file lo-math.h
John W. Eaton <jwe@octave.org>
parents:
23649
diff
changeset
|
88 |
5742 | 89 #include "lo-ieee.h" |
23475
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
90 #include "randgamma.h" |
5742 | 91 #include "randmtzig.h" |
92 | |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
93 namespace octave |
5742 | 94 { |
25434
859ad1f0b85e
use templates to eliminate duplicate code in random number generators
John W. Eaton <jwe@octave.org>
parents:
25433
diff
changeset
|
95 template <typename T> void rand_gamma (T a, octave_idx_type n, T *r) |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
96 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
97 octave_idx_type i; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
98 /* If a < 1, start by generating gamma (1+a) */ |
25434
859ad1f0b85e
use templates to eliminate duplicate code in random number generators
John W. Eaton <jwe@octave.org>
parents:
25433
diff
changeset
|
99 const T d = (a < 1. ? 1.+a : a) - 1./3.; |
859ad1f0b85e
use templates to eliminate duplicate code in random number generators
John W. Eaton <jwe@octave.org>
parents:
25433
diff
changeset
|
100 const T c = 1./std::sqrt (9.*d); |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
101 |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
102 /* Handle invalid cases */ |
25434
859ad1f0b85e
use templates to eliminate duplicate code in random number generators
John W. Eaton <jwe@octave.org>
parents:
25433
diff
changeset
|
103 if (a <= 0 || lo_ieee_isinf (a)) |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
104 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
105 for (i=0; i < n; i++) |
25434
859ad1f0b85e
use templates to eliminate duplicate code in random number generators
John W. Eaton <jwe@octave.org>
parents:
25433
diff
changeset
|
106 r[i] = numeric_limits<T>::NaN (); |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
107 return; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
108 } |
5742 | 109 |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
110 for (i=0; i < n; i++) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
111 { |
25434
859ad1f0b85e
use templates to eliminate duplicate code in random number generators
John W. Eaton <jwe@octave.org>
parents:
25433
diff
changeset
|
112 T x, xsq, v, u; |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
113 restart: |
25434
859ad1f0b85e
use templates to eliminate duplicate code in random number generators
John W. Eaton <jwe@octave.org>
parents:
25433
diff
changeset
|
114 x = rand_normal<T> (); |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
115 v = (1+c*x); |
26437
b43eb366666c
randgamma.cc: Fix static analyzer detected issues (bug #55347).
Rik <rik@octave.org>
parents:
26376
diff
changeset
|
116 v *= (v*v); |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
117 if (v <= 0) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
118 goto restart; /* rare, so don't bother moving up */ |
25434
859ad1f0b85e
use templates to eliminate duplicate code in random number generators
John W. Eaton <jwe@octave.org>
parents:
25433
diff
changeset
|
119 u = rand_uniform<T> (); |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
120 xsq = x*x; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
121 if (u >= 1.-0.0331*xsq*xsq && std::log (u) >= 0.5*xsq + d*(1-v+std::log (v))) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
122 goto restart; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
123 r[i] = d*v; |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
124 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
125 if (a < 1) |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
126 { |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
127 /* Use gamma(a) = gamma(1+a)*U^(1/a) */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
128 /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */ |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
129 for (i = 0; i < n; i++) |
25434
859ad1f0b85e
use templates to eliminate duplicate code in random number generators
John W. Eaton <jwe@octave.org>
parents:
25433
diff
changeset
|
130 r[i] *= exp (-rand_exponential<T> () / a); |
25433
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
131 } |
49e0447413ad
use templates and move rand functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
25054
diff
changeset
|
132 } |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
133 |
29228
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27923
diff
changeset
|
134 template OCTAVE_API void rand_gamma (double, octave_idx_type, double *); |
5c14f81e0937
Set API tags in files in liboctave/numeric (patch #8919).
Markus Mützel <markus.muetzel@gmx.de>
parents:
27923
diff
changeset
|
135 template OCTAVE_API void rand_gamma (float, octave_idx_type, float *); |
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
136 } |