diff options
Diffstat (limited to 'numpy/random/mtrand/distributions.c')
-rw-r--r-- | numpy/random/mtrand/distributions.c | 82 |
1 files changed, 41 insertions, 41 deletions
diff --git a/numpy/random/mtrand/distributions.c b/numpy/random/mtrand/distributions.c index d792cf86d..05f62f0cf 100644 --- a/numpy/random/mtrand/distributions.c +++ b/numpy/random/mtrand/distributions.c @@ -7,10 +7,10 @@ * distribute, sublicense, and/or sell copies of the Software, and to * permit persons to whom the Software is furnished to do so, subject to * the following conditions: - * + * * The above copyright notice and this permission notice shall be included * in all copies or substantial portions of the Software. - * + * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. @@ -21,9 +21,9 @@ */ /* The implementations of rk_hypergeometric_hyp(), rk_hypergeometric_hrua(), - * and rk_triangular() were adapted from Ivan Frohne's rv.py which has this + * and rk_triangular() were adapted from Ivan Frohne's rv.py which has this * license: - * + * * Copyright 1998 by Ivan Frohne; Wasilla, Alaska, U.S.A. * All Rights Reserved * @@ -52,8 +52,8 @@ #ifndef M_PI #define M_PI 3.14159265358979323846264338328 -#endif -/* log-gamma function to support some of these distributions. The +#endif +/* log-gamma function to support some of these distributions. The * algorithm comes from SPECFUN by Shanjie Zhang and Jianming Jin and their * book "Computation of Special Functions", 1996, John Wiley & Sons, Inc. */ @@ -62,7 +62,7 @@ double loggam(double x) { double x0, x2, xp, gl, gl0; long k, n; - + static double a[10] = {8.333333333333333e-02,-2.777777777777778e-03, 7.936507936507937e-04,-5.952380952380952e-04, 8.417508417508418e-04,-1.917526917526918e-03, @@ -162,7 +162,7 @@ double rk_standard_gamma(rk_state *state, double shape) { do { - X = rk_gauss(state); + X = rk_gauss(state); V = 1.0 + c*X; } while (V <= 0.0); @@ -225,7 +225,7 @@ double rk_noncentral_chisquare(rk_state *state, double df, double nonc) double rk_f(rk_state *state, double dfnum, double dfden) { - return ((rk_chisquare(state, dfnum) * dfden) / + return ((rk_chisquare(state, dfnum) * dfden) / (rk_chisquare(state, dfden) * dfnum)); } @@ -241,7 +241,7 @@ long rk_binomial_btpe(rk_state *state, long n, double p) double a,u,v,s,F,rho,t,A,nrq,x1,x2,f1,f2,z,z2,w,w2,x; long m,y,k,i; - if (!(state->has_binomial) || + if (!(state->has_binomial) || (state->nsave != n) || (state->psave != p)) { @@ -380,7 +380,7 @@ long rk_binomial_inversion(rk_state *state, long n, double p) double q, qn, np, px, U; long X, bound; - if (!(state->has_binomial) || + if (!(state->has_binomial) || (state->nsave != n) || (state->psave != p)) { @@ -404,12 +404,12 @@ long rk_binomial_inversion(rk_state *state, long n, double p) while (U > px) { X++; - if (X > bound) + if (X > bound) { X = 0; px = qn; U = rk_double(state); - } else + } else { U -= px; px = ((n-X+1) * p * px)/(X*q); @@ -514,7 +514,7 @@ long rk_poisson_ptrs(rk_state *state, double lam) return k; } - + } } @@ -525,11 +525,11 @@ long rk_poisson(rk_state *state, double lam) { return rk_poisson_ptrs(state, lam); } - else if (lam == 0) + else if (lam == 0) { return 0; } - else + else { return rk_poisson_mult(state, lam); } @@ -551,7 +551,7 @@ double rk_standard_t(rk_state *state, double df) } /* Uses the rejection algorithm compared against the wrapped Cauchy - distribution suggested by Best and Fisher and documented in + distribution suggested by Best and Fisher and documented in Chapter 9 of Luc's Non-Uniform Random Variate Generation. http://cg.scs.carleton.ca/~luc/rnbookindex.html (but corrected to match the algorithm in R and Python) @@ -600,7 +600,7 @@ double rk_vonmises(rk_state *state, double mu, double kappa) if (neg) { mod *= -1; - } + } return mod; } @@ -624,12 +624,12 @@ double rk_power(rk_state *state, double a) double rk_laplace(rk_state *state, double loc, double scale) { double U; - + U = rk_double(state); if (U < 0.5) { U = loc + scale * log(U + U); - } else + } else { U = loc - scale * log(2.0 - U - U); } @@ -639,7 +639,7 @@ double rk_laplace(rk_state *state, double loc, double scale) double rk_gumbel(rk_state *state, double loc, double scale) { double U; - + U = 1.0 - rk_double(state); return loc - scale * log(-log(U)); } @@ -647,7 +647,7 @@ double rk_gumbel(rk_state *state, double loc, double scale) double rk_logistic(rk_state *state, double loc, double scale) { double U; - + U = rk_double(state); return loc + scale * log(U/(1.0 - U)); } @@ -666,7 +666,7 @@ double rk_wald(rk_state *state, double mean, double scale) { double U, X, Y; double mu_2l; - + mu_2l = mean / (2*scale); Y = rk_gauss(state); Y = mean*Y*Y; @@ -710,7 +710,7 @@ long rk_geometric_search(rk_state *state, double p) double U; long X; double sum, prod, q; - + X = 1; sum = prod = p; q = 1.0 - p; @@ -744,10 +744,10 @@ long rk_hypergeometric_hyp(rk_state *state, long good, long bad, long sample) { long d1, K, Z; double d2, U, Y; - + d1 = bad + good - sample; d2 = (double)min(bad, good); - + Y = d2; K = sample; while (Y > 0.0) @@ -772,7 +772,7 @@ long rk_hypergeometric_hrua(rk_state *state, long good, long bad, long sample) double d4, d5, d6, d7, d8, d10, d11; long Z; double T, W, X, Y; - + mingoodbad = min(good, bad); popsize = good + bad; maxgoodbad = max(good, bad); @@ -783,39 +783,39 @@ long rk_hypergeometric_hrua(rk_state *state, long good, long bad, long sample) d7 = sqrt((popsize - m) * sample * d4 *d5 / (popsize-1) + 0.5); d8 = D1*d7 + D2; d9 = (long)floor((double)((m+1)*(mingoodbad+1))/(popsize+2)); - d10 = (loggam(d9+1) + loggam(mingoodbad-d9+1) + loggam(m-d9+1) + + d10 = (loggam(d9+1) + loggam(mingoodbad-d9+1) + loggam(m-d9+1) + loggam(maxgoodbad-m+d9+1)); d11 = min(min(m, mingoodbad)+1.0, floor(d6+16*d7)); /* 16 for 16-decimal-digit precision in D1 and D2 */ - + while (1) { X = rk_double(state); Y = rk_double(state); W = d6 + d8*(Y- 0.5)/X; - + /* fast rejection: */ if ((W < 0.0) || (W >= d11)) continue; - + Z = (long)floor(W); T = d10 - (loggam(Z+1) + loggam(mingoodbad-Z+1) + loggam(m-Z+1) + loggam(maxgoodbad-m+Z+1)); - + /* fast acceptance: */ if ((X*(4.0-X)-3.0) <= T) break; - + /* fast rejection: */ if (X*(X-T) >= 1) continue; if (2.0*log(X) <= T) break; /* acceptance */ } - + /* this is a correction to HRUA* by Ivan Frohne in rv.py */ if (good > bad) Z = m - Z; - + /* another fix from rv.py to allow sample to exceed popsize/2 */ if (m < sample) Z = good - Z; - + return Z; } #undef D1 @@ -836,20 +836,20 @@ double rk_triangular(rk_state *state, double left, double mode, double right) { double base, leftbase, ratio, leftprod, rightprod; double U; - + base = right - left; leftbase = mode - left; ratio = leftbase / base; leftprod = leftbase*base; rightprod = (right - mode)*base; - + U = rk_double(state); if (U <= ratio) { return left + sqrt(U*leftprod); - } else + } else { - return right - sqrt((1.0 - U) * rightprod); + return right - sqrt((1.0 - U) * rightprod); } } @@ -857,7 +857,7 @@ long rk_logseries(rk_state *state, double p) { double q, r, U, V; long result; - + r = log(1.0 - p); while (1) { |