diff options
Diffstat (limited to 'numpy/random/mtrand/distributions.c')
| -rw-r--r-- | numpy/random/mtrand/distributions.c | 845 |
1 files changed, 845 insertions, 0 deletions
diff --git a/numpy/random/mtrand/distributions.c b/numpy/random/mtrand/distributions.c new file mode 100644 index 000000000..3f5ff2355 --- /dev/null +++ b/numpy/random/mtrand/distributions.c @@ -0,0 +1,845 @@ +/* Copyright 2005 Robert Kern (robert.kern@gmail.com) + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * 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. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +/* The implementations of rk_hypergeometric_hyp(), rk_hypergeometric_hrua(), + * 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 + * + * Permission to use, copy, modify and distribute this software and its + * documentation for any purpose, free of charge, is granted 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 AND DOCUMENTATION IS PROVIDED WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO MERCHANTABILITY, FITNESS + * FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR + * OR COPYRIGHT HOLDER BE LIABLE FOR ANY CLAIM OR DAMAGES IN A CONTRACT + * ACTION, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR ITS DOCUMENTATION. + */ + +#include <math.h> +#include "distributions.h" +#include <stdio.h> + +#ifndef min +#define min(x,y) ((x<y)?x:y) +#define max(x,y) ((x>y)?x:y) +#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. + */ +extern double loggam(double x); +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, + 6.410256410256410e-03,-2.955065359477124e-02, + 1.796443723688307e-01,-1.39243221690590e+00}; + x0 = x; + n = 0; + if ((x == 1.0) || (x == 2.0)) + { + return 0.0; + } + else if (x <= 7.0) + { + n = (long)(7 - x); + x0 = x + n; + } + x2 = 1.0/(x0*x0); + xp = 2*M_PI; + gl0 = a[9]; + for (k=8; k>=0; k--) + { + gl0 *= x2; + gl0 += a[k]; + } + gl = gl0/x0 + 0.5*log(xp) + (x0-0.5)*log(x0) - x0; + if (x <= 7.0) + { + for (k=1; k<=n; k++) + { + gl -= log(x0-1.0); + x0 -= 1.0; + } + } + return gl; +} + +double rk_normal(rk_state *state, double loc, double scale) +{ + return loc + scale*rk_gauss(state); +} + +double rk_standard_exponential(rk_state *state) +{ + /* We use -log(1-U) since U is [0, 1) */ + return -log(1.0 - rk_double(state)); +} + +double rk_exponential(rk_state *state, double scale) +{ + return scale * rk_standard_exponential(state); +} + +double rk_uniform(rk_state *state, double loc, double scale) +{ + return loc + scale*rk_double(state); +} + +double rk_standard_gamma(rk_state *state, double shape) +{ + double b, c; + double U, V, X, Y; + + if (shape == 1.0) + { + return rk_standard_exponential(state); + } + else if (shape < 1.0) + { + for (;;) + { + U = rk_double(state); + V = rk_standard_exponential(state); + if (U <= 1.0 - shape) + { + X = pow(U, 1./shape); + if (X <= V) + { + return X; + } + } + else + { + Y = -log((1-U)/shape); + X = pow(1.0 - shape + shape*Y, 1./shape); + if (X <= (V + Y)) + { + return X; + } + } + } + } + else + { + b = shape - 1./3.; + c = 1./sqrt(9*b); + for (;;) + { + do + { + X = rk_gauss(state); + V = 1.0 + c*X; + } while (V <= 0.0); + + V = V*V*V; + U = rk_double(state); + if (U < 1.0 - 0.0331*(X*X)*(X*X)) return (b*V); + if (log(U) < 0.5*X*X + b*(1. - V + log(V))) return (b*V); + } + } +} + +double rk_gamma(rk_state *state, double shape, double scale) +{ + return scale * rk_standard_gamma(state, shape); +} + +double rk_beta(rk_state *state, double a, double b) +{ + double Ga, Gb; + + if ((a <= 1.0) && (b <= 1.0)) + { + double U, V, X, Y; + /* Use Jonk's algorithm */ + + while (1) + { + U = rk_double(state); + V = rk_double(state); + X = pow(U, 1.0/a); + Y = pow(V, 1.0/b); + + if ((X + Y) <= 1.0) + { + return X; + } + } + } + else + { + Ga = rk_standard_gamma(state, a); + Gb = rk_standard_gamma(state, b); + return Ga/(Ga + Gb); + } +} + +double rk_chisquare(rk_state *state, double df) +{ + return 2.0*rk_standard_gamma(state, df/2.0); +} + +double rk_noncentral_chisquare(rk_state *state, double df, double nonc) +{ + double Chi2, N; + + Chi2 = rk_chisquare(state, df-1); + N = rk_gauss(state) + sqrt(nonc); + return Chi2 + N*N; +} + +double rk_f(rk_state *state, double dfnum, double dfden) +{ + return rk_chisquare(state, dfnum) / rk_chisquare(state, dfden); +} + +double rk_noncentral_f(rk_state *state, double dfnum, double dfden, double nonc) +{ + return ((rk_noncentral_chisquare(state, dfnum, nonc)*dfden) / + (rk_chisquare(state, dfden)*dfnum)); +} + +long rk_binomial_btpe(rk_state *state, long n, double p) +{ + double r,q,fm,p1,xm,xl,xr,c,laml,lamr,p2,p3,p4; + 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) || + (state->nsave != n) || + (state->psave != p)) + { + /* initialize */ + state->nsave = n; + state->psave = p; + state->has_binomial = 1; + state->r = r = min(p, 1.0-p); + state->q = q = 1.0 - r; + state->fm = fm = n*r+r; + state->m = m = (long)floor(state->fm); + state->p1 = p1 = floor(2.195*sqrt(n*r*q)-4.6*q) + 0.5; + state->xm = xm = m + 0.5; + state->xl = xl = xm - p1; + state->xr = xr = xm + p1; + state->c = c = 0.134 + 20.5/(15.3 + m); + a = (fm - xl)/(fm-xl*r); + state->laml = laml = a*(1.0 + a/2.0); + a = (xr - fm)/(xr*q); + state->lamr = lamr = a*(1.0 + a/2.0); + state->p2 = p2 = p1*(1.0 + 2.0*c); + state->p3 = p3 = p2 + c/laml; + state->p4 = p4 = p3 + c/lamr; + } + else + { + r = state->r; + q = state->q; + fm = state->fm; + m = state->m; + p1 = state->p1; + xm = state->xm; + xl = state->xl; + xr = state->xr; + c = state->c; + laml = state->laml; + lamr = state->lamr; + p2 = state->p2; + p3 = state->p3; + p4 = state->p4; + } + + /* sigh ... */ + Step10: + nrq = n*r*q; + u = rk_double(state)*p4; + v = rk_double(state); + if (u > p1) goto Step20; + y = (long)floor(xm - p1*v + u); + goto Step60; + + Step20: + if (u > p2) goto Step30; + x = xl + (u - p1)/c; + v = v*c + 1.0 - fabs(m - x + 0.5)/p1; + if (v > 1.0) goto Step10; + y = (long)floor(x); + goto Step50; + + Step30: + if (u > p3) goto Step40; + y = (long)floor(xl + log(v)/laml); + if (y < 0) goto Step10; + v = v*(u-p2)*laml; + goto Step50; + + Step40: + y = (int)floor(xr - log(v)/lamr); + if (y > n) goto Step10; + v = v*(u-p3)*lamr; + + Step50: + k = fabs(y - m); + if ((k > 20) && (k < ((nrq)/2.0 - 1))) goto Step52; + + s = r/q; + a = s*(n+1); + F = 1.0; + if (m < y) + { + for (i=m; i<=y; i++) + { + F *= (a/i - s); + } + } + else if (m > y) + { + for (i=y; i<=m; i++) + { + F /= (a/i - s); + } + } + else + { + if (v > F) goto Step10; + goto Step60; + } + + Step52: + rho = (k/(nrq))*((k*(k/3.0 + 0.625) + 0.16666666666666666)/nrq + 0.5); + t = -k*k/(2*nrq); + A = log(v); + if (A < (t - rho)) goto Step60; + if (A > (t + rho)) goto Step10; + + x1 = y+1; + f1 = m+1; + z = n+1-m; + w = n-y+1; + x2 = x1*x1; + f2 = f1*f1; + z2 = z*z; + w2 = w*w; + if (A > (xm*log(f1/x1) + + (n-m+0.5)*log(z/w) + + (y-m)*log(w*r/(x1*q)) + + (13680.-(462.-(132.-(99.-140./f2)/f2)/f2)/f2)/f1/166320. + + (13680.-(462.-(132.-(99.-140./z2)/z2)/z2)/z2)/z/166320. + + (13680.-(462.-(132.-(99.-140./x2)/x2)/x2)/x2)/x1/166320. + + (13680.-(462.-(132.-(99.-140./w2)/w2)/w2)/w2)/w/166320.)) + { + goto Step10; + } + + Step60: + if (p > 0.5) + { + y = n - y; + } + + return y; +} + +long rk_binomial_inversion(rk_state *state, long n, double p) +{ + double q, qn, np, px, U; + long X, bound; + + if (!(state->has_binomial) || + (state->nsave != n) || + (state->psave != p)) + { + state->nsave = n; + state->psave = p; + state->has_binomial = 1; + state->q = q = 1.0 - p; + state->r = qn = exp(n * log(q)); + state->c = np = n*p; + state->m = bound = min(n, np + 10.0*sqrt(np)); + } else + { + q = state->q; + qn = state->r; + np = state->c; + bound = state->m; + } + X = 0; + px = qn; + U = rk_double(state); + while (U > px) + { + X++; + if (X > bound) + { + X = 0; + px = qn; + U = rk_double(state); + } else + { + U -= px; + px = ((n-X+1) * p * px)/(X*q); + } + } + return X; +} + +long rk_binomial(rk_state *state, long n, double p) +{ + double q; + + if (p <= 0.5) + { + if (p*n <= 30.0) + { + return rk_binomial_inversion(state, n, p); + } + else + { + return rk_binomial_btpe(state, n, p); + } + } + else + { + q = 1.0-p; + if (q*n <= 30.0) + { + return n - rk_binomial_inversion(state, n, q); + } + else + { + return n - rk_binomial_btpe(state, n, q); + } + } + +} + +long rk_negative_binomial(rk_state *state, long n, double p) +{ + double Y; + + Y = rk_gamma(state, n, (1-p)/p); + return rk_poisson(state, Y); +} + +long rk_poisson_mult(rk_state *state, double lam) +{ + long X; + double prod, U, enlam; + + enlam = exp(-lam); + X = 0; + prod = 1.0; + while (1) + { + U = rk_double(state); + prod *= U; + if (prod > enlam) + { + X += 1; + } + else + { + return X; + } + } +} + +#define LS2PI 0.91893853320467267 +#define TWELFTH 0.083333333333333333333333 +long rk_poisson_ptrs(rk_state *state, double lam) +{ + long k; + double U, V, slam, loglam, a, b, invalpha, vr, us; + + slam = sqrt(lam); + loglam = log(lam); + b = 0.931 + 2.53*slam; + a = -0.059 + 0.02483*b; + invalpha = 1.1239 + 1.1328/(b-3.4); + vr = 0.9277 - 3.6224/(b-2); + + while (1) + { + U = rk_double(state) - 0.5; + V = rk_double(state); + us = 0.5 - fabs(U); + k = (long)floor((2*a/us + b)*U + lam + 0.43); + if ((us >= 0.07) && (V <= vr)) + { + return k; + } + if ((k < 0) || + ((us < 0.013) && (V > us))) + { + continue; + } + if ((log(V) + log(invalpha) - log(a/(us*us)+b)) <= + (-lam + k*loglam - loggam(k+1))) + { + return k; + } + + + } + +} + +long rk_poisson(rk_state *state, double lam) +{ + if (lam >= 10) + { + return rk_poisson_ptrs(state, lam); + } + else + { + return rk_poisson_mult(state, lam); + } +} + +double rk_standard_cauchy(rk_state *state) +{ + return rk_gauss(state) / rk_gauss(state); +} + +double rk_standard_t(rk_state *state, double df) +{ + double N, G, X; + + N = rk_gauss(state); + G = rk_standard_gamma(state, df/2); + X = sqrt(df/2)*N/sqrt(G); + return X; +} + +double rk_vonmises(rk_state *state, double mu, double kappa) +{ + double r, rho, s; + double U, V, W, Y, Z; + double result, mod; + + if (kappa < 1e-8) + { + return M_PI * (2*rk_double(state)-1); + } + else + { + r = 1 + sqrt(1 + 4*kappa*kappa); + rho = (r - sqrt(2*r))/(2*kappa); + s = (1 + rho*rho)/(2*rho); + + while (1) + { + U = 2*rk_double(state) - 1; + V = 2*rk_double(state) - 1; + Z = cos(M_PI*U); + W = (1 + s*Z)/(s + Z); + Y = kappa * (s - W); + if ((Y*(2-Y) - V >= 0) || (log(Y/V)+1 - Y >= 0)) + { + break; + } + } + + if (U < 0) + { + result = acos(W); + } + else + { + result = -acos(W); + } + result += mu + M_PI; + mod = fmod(result, 2*M_PI); + if (mod && (mod < 0)) + { + mod += 2*M_PI; + } + return mod - M_PI; + } +} + +double rk_pareto(rk_state *state, double a) +{ + return exp(rk_standard_exponential(state)/a) - 1; +} + +double rk_weibull(rk_state *state, double a) +{ + return pow(rk_standard_exponential(state), 1./a); +} + +double rk_power(rk_state *state, double a) +{ + return pow(1 - exp(-rk_standard_exponential(state)), 1./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 + { + U = loc - scale * log(2.0 - U - U); + } + return U; +} + +double rk_gumbel(rk_state *state, double loc, double scale) +{ + double U; + + U = 1.0 - rk_double(state); + return loc - scale * log(-log(U)); +} + +double rk_logistic(rk_state *state, double loc, double scale) +{ + double U; + + U = rk_double(state); + return loc + scale * log(U/(1.0 - U)); +} + +double rk_lognormal(rk_state *state, double mean, double sigma) +{ + return exp(rk_normal(state, mean, sigma)); +} + +double rk_rayleigh(rk_state *state, double mode) +{ + return mode*sqrt(-2.0 * log(1.0 - rk_double(state))); +} + +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; + X = mean + mu_2l*(Y - sqrt(4*scale*Y + Y*Y)); + U = rk_double(state); + if (U <= mean/(mean+X)) + { + return X; + } else + { + return mean*mean/X; + } +} + +long rk_zipf(rk_state *state, double a) +{ + double T, U, V; + long X; + double b; + + b = pow(2.0, a-1.0); + do + { + U = rk_double(state); + V = rk_double(state); + X = (long)floor(pow(U, -1.0/(a-1.0))); + T = pow(1.0 + 1.0/X, a-1.0); + } while ((V *X*(T-1.0)/(b-1.0)) > (T/b)); + return X; +} + +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; + U = rk_double(state); + while (U > sum) + { + prod *= q; + sum += prod; + X++; + } + return X; +} + +long rk_geometric_inversion(rk_state *state, double p) +{ + return (long)ceil(log(1.0-rk_double(state))/log(1.0-p)); +} + +long rk_geometric(rk_state *state, double p) +{ + if (p >= 0.333333333333333333333333) + { + return rk_geometric_search(state, p); + } else + { + return rk_geometric_inversion(state, p); + } +} + +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) + { + U = rk_double(state); + Y -= (long)floor(U + Y/(d1 + K)); + K--; + if (K == 0) break; + } + Z = (long)(d2 - Y); + if (bad > good) Z = sample - Z; + return Z; +} + +/* D1 = 2*sqrt(2/e) */ +/* D2 = 3 - 2*sqrt(3/e) */ +#define D1 1.7155277699214135 +#define D2 0.8989161620588988 +long rk_hypergeometric_hrua(rk_state *state, long good, long bad, long sample) +{ + long mingoodbad, maxgoodbad, popsize, m, d9; + 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); + m = min(sample, popsize - sample); + d4 = ((double)mingoodbad) / popsize; + d5 = 1.0 - d4; + d6 = m*d4 + 0.5; + 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) + + 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 (bad > good) Z = m - Z; + + /* another fix from rv.py to allow sample to exceed popsize/2 */ + if (m < sample) Z = bad - Z; + + return Z; +} +#undef D1 +#undef D2 + +long rk_hypergeometric(rk_state *state, long good, long bad, long sample) +{ + if (sample > 10) + { + return rk_hypergeometric_hrua(state, good, bad, sample); + } else + { + return rk_hypergeometric_hyp(state, good, bad, sample); + } +} + +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 + { + return right - sqrt((1.0 - U) * rightprod); + } +} + +long rk_logseries(rk_state *state, double p) +{ + double q, r, U, V; + + r = log(1.0 - p); + + V = rk_double(state); + if (V >= p) return 1; + U = rk_double(state); + q = 1.0 - exp(r*U); + if (V <= q*q) return (long)floor(1 + log(V)/log(q)); + if (V <= q) return 1; + return 2; +} |