octave: liboctave/numeric/randgamma.cc annotate

annotate liboctave/numeric/randgamma.cc @ 31249:de6fc38c78c6

Make Jacobian types offered by dlsode.f accessible by lsode (bug #31626). * liboctave/numeric/LSODE-opts.in: Add options "jacobian type", "lower jacobian subdiagonals", and "upper jacobian subdiagonals". * liboctave/numeric/LSODE.cc (file scope, lsode_j, LSODE::do_integrate (double)): Handle new configurable Jacobian types. * build-aux/mk-opts.pl: Don't implicitly convert to integer in condition.

author	Olaf Till <olaf.till@uni-jena.de>
date	Fri, 12 Nov 2010 08:53:05 +0100
parents	796f54d4ddbf
children	e88a07dec498

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 //
30564 796f54d4ddbf update Octave Project Developers copyright for the new year John W. Eaton <jwe@octave.org> parents: 29359 diff changeset	3 // Copyright (C) 2006-2022 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 4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	25
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	26 /* Original version written by Paul Kienzle distributed as free
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe] jwe parents: 5742 diff changeset	27 software in the in the public domain. */
5742 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	28
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	29 /*
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	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 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	33
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	34 Generate a series of standard gamma distributions.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	35
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	36 See: Marsaglia G and Tsang W (2000), "A simple method for generating
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	37 gamma variables", ACM Transactions on Mathematical Software 26(3) 363-372
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	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 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	40 * NAN: value to return for Not-A-Number
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	41 * RUNI: uniform generator on (0,1)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	42 * RNOR: normal generator
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	43 * REXP: exponential generator, or -log(RUNI) if one isn't available
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	44 * INFINITE: function to test whether a value is infinite
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	45
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	46 Test using:
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	47 mean = a
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	48 variance = a
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	49 skewness = 2/sqrt(a)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	50 kurtosis = 3 + 6/sqrt(a)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	51
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	52 Note that randg can be used to generate many distributions:
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	53
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	54 gamma(a,b) for a>0, b>0 (from R)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	55 r = b*randg(a)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	56 beta(a,b) for a>0, b>0
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	57 r1 = randg(a,1)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	58 r = r1 / (r1 + randg(b,1))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	59 Erlang(a,n)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	60 r = a*randg(n)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	61 chisq(df) for df>0
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	62 r = 2*randg(df/2)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	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 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	65 F(n1,n2) for 0<n1, 0<n2
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	66 r1 = 2*randg(n1/2)/n1 or 1 if n1 is infinite
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	67 r2 = 2*randg(n2/2)/n2 or 1 if n2 is infinite
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	68 r = r1 / r2
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	69 negative binonial (n, p) for n>0, 0<p<=1
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	70 r = randp((1-p)/p * randg(n))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	71 (from R, citing Devroye(1986), Non-Uniform Random Variate Generation)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	72 non-central chisq(df,L), for df>=0 and L>0 (use chisq if L=0)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	73 r = randp(L/2)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	74 r(r>0) = 2*randg(r(r>0))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	75 r(df>0) += 2*randg(df(df>0)/2)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	76 (from R, citing formula 29.5b-c in Johnson, Kotz, Balkrishnan(1995))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	77 Dirichlet(a1,...,ak) for ai > 0
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	78 r = (randg(a1),...,randg(ak))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	79 r = r / sum(r)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	80 (from GSL, citing Law & Kelton(1991), Simulation Modeling and Analysis)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	81 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	82
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	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 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	85 #endif
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	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 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	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 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	91 #include "randmtzig.h"
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	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 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	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 2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe] jwe parents: diff changeset	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 = (vv);
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.0331xsqxsq && std::log (u) >= 0.5xsq + 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 }

Mercurial > octave

annotate liboctave/numeric/randgamma.cc @ 31249:de6fc38c78c6