The R Project SVN R

Rev

Rev 70052 | Rev 77694 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
574 ihaka 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
3 
 *  Copyright (C) 1998 Ross Ihaka
68749 maechler 4 
 *  Copyright (C) 2000-2015 The R Core Team
5 
 *  Copyright (C) 2005-2015 The R Foundation
574 ihaka 6 
 *
7 
 *  This program is free software; you can redistribute it and/or modify
8 
 *  it under the terms of the GNU General Public License as published by
9 
 *  the Free Software Foundation; either version 2 of the License, or
10 
 *  (at your option) any later version.
11 
 *
12 
 *  This program is distributed in the hope that it will be useful,
13 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 
 *  GNU General Public License for more details.
16 
 *
42302 ripley 17 
 *  You should have received a copy of the GNU General Public License
42310 ripley 18 
 *  along with this program; if not, a copy is available at
68947 ripley 19 
 *  https://www.R-project.org/Licenses/
574 ihaka 20 
 *
21 
 *  SYNOPSIS
22 
 *
19500 hornik 23 
 *    #include <Rmath.h>
574 ihaka 24 
 *    double rhyper(double NR, double NB, double n);
25 
 *
26 
 *  DESCRIPTION
27 
 *
28 
 *    Random variates from the hypergeometric distribution.
29 
 *    Returns the number of white balls drawn when kk balls
30 
 *    are drawn at random from an urn containing nn1 white
31 
 *    and nn2 black balls.
32 
 *
33 
 *  REFERENCE
34 
 *
35 
 *    V. Kachitvichyanukul and B. Schmeiser (1985).
36 
 *    ``Computer generation of hypergeometric random variates,''
37 
 *    Journal of Statistical Computation and Simulation 22, 127-145.
33718 maechler 38 
 *
39 
 *    The original algorithm had a bug -- R bug report PR#7314 --
40 
 *    giving numbers slightly too small in case III h2pe
41 
 *    where (m < 100 || ix <= 50) , see below.
574 ihaka 42 
 */
43 
 
8431 ripley 44 
#include "nmath.h"
64655 ripley 45 
#include "dpq.h"
70052 maechler 46 
#include <limits.h>
574 ihaka 47 
 
70052 maechler 48 
// afc(i) :=  ln( i! )	[logarithm of the factorial i]
574 ihaka 49 
static double afc(int i)
50 
{
70052 maechler 51 
    // If (i > 7), use Stirling's approximation, otherwise use table lookup.
52 
    const static double al[8] =
13434 maechler 53 
    {
54 
	0.0,/*ln(0!)=ln(1)*/
55 
	0.0,/*ln(1!)=ln(1)*/
56 
	0.69314718055994530941723212145817,/*ln(2) */
57 
	1.79175946922805500081247735838070,/*ln(6) */
58 
	3.17805383034794561964694160129705,/*ln(24)*/
59 
	4.78749174278204599424770093452324,
60 
	6.57925121201010099506017829290394,
61 
	8.52516136106541430016553103634712
70052 maechler 62 
	/* 10.60460290274525022841722740072165, approx. value below =
63 
	   10.6046028788027; rel.error = 2.26 10^{-9}
64 
 
65 
	  FIXME: Use constants and if(n > ..) decisions from ./stirlerr.c
66 
	  -----  will be even *faster* for n > 500 (or so)
67 
	*/
13434 maechler 68 
    };
69 
 
1028 maechler 70 
    if (i < 0) {
70052 maechler 71 
	MATHLIB_WARNING(("rhyper.c: afc(i), i=%d < 0 -- SHOULD NOT HAPPEN!\n"), i);
72 
	return -1; // unreached
574 ihaka 73 
    }
70052 maechler 74 
    if (i <= 7)
75 
	return al[i];
76 
    // else i >= 8 :
77 
    double di = i, i2 = di*di;
78 
    return (di + 0.5) * log(di) - di + M_LN_SQRT_2PI +
79 
	(0.0833333333333333 - 0.00277777777777778 / i2) / di;
574 ihaka 80 
}
81 
 
70052 maechler 82 
//     rhyper(NR, NB, n) -- NR 'red', NB 'blue', n drawn, how many are 'red'
574 ihaka 83 
double rhyper(double nn1in, double nn2in, double kkin)
84 
{
13434 maechler 85 
    /* extern double afc(int); */
86 
 
574 ihaka 87 
    int nn1, nn2, kk;
70052 maechler 88 
    int ix; // return value (coerced to double at the very end)
89 
    Rboolean setup1, setup2;
33718 maechler 90 
 
68749 maechler 91 
    /* These should become 'thread_local globals' : */
92 
    static int ks = -1, n1s = -1, n2s = -1;
70052 maechler 93 
    static int m, minjx, maxjx;
94 
    static int k, n1, n2; // <- not allowing larger integer par
68757 maechler 95 
    static double tn;
70052 maechler 96 
 
68749 maechler 97 
    // II :
98 
    static double w;
99 
    // III:
100 
    static double a, d, s, xl, xr, kl, kr, lamdl, lamdr, p1, p2, p3;
574 ihaka 101 
 
102 
    /* check parameter validity */
103 
 
7671 maechler 104 
    if(!R_FINITE(nn1in) || !R_FINITE(nn2in) || !R_FINITE(kkin))
105 
	ML_ERR_return_NAN;
574 ihaka 106 
 
70052 maechler 107 
    nn1in = R_forceint(nn1in);
108 
    nn2in = R_forceint(nn2in);
109 
    kkin  = R_forceint(kkin);
574 ihaka 110 
 
70052 maechler 111 
    if (nn1in < 0 || nn2in < 0 || kkin < 0 || kkin > nn1in + nn2in)
7671 maechler 112 
	ML_ERR_return_NAN;
70052 maechler 113 
    if (nn1in >= INT_MAX || nn2in >= INT_MAX || kkin >= INT_MAX) {
114 
	/* large n -- evade integer overflow (and inappropriate algorithms)
115 
	   -------- */
116 
        // FIXME: Much faster to give rbinom() approx when appropriate; -> see Kuensch(1989)
117 
	// Johnson, Kotz,.. p.258 (top) mention the *four* different binomial approximations
118 
	if(kkin == 1.) { // Bernoulli
119 
	    return rbinom(kkin, nn1in / (nn1in + nn2in));
120 
	}
121 
	// Slow, but safe: return  F^{-1}(U)  where F(.) = phyper(.) and  U ~ U[0,1]
122 
	return qhyper(unif_rand(), nn1in, nn2in, kkin, FALSE, FALSE);
123 
    }
124 
    nn1 = (int)nn1in;
125 
    nn2 = (int)nn2in;
126 
    kk  = (int)kkin;
7671 maechler 127 
 
574 ihaka 128 
    /* if new parameter values, initialize */
129 
    if (nn1 != n1s || nn2 != n2s) {
13434 maechler 130 
	setup1 = TRUE;	setup2 = TRUE;
574 ihaka 131 
    } else if (kk != ks) {
13434 maechler 132 
	setup1 = FALSE;	setup2 = TRUE;
133 
    } else {
134 
	setup1 = FALSE;	setup2 = FALSE;
574 ihaka 135 
    }
136 
    if (setup1) {
137 
	n1s = nn1;
138 
	n2s = nn2;
139 
	tn = nn1 + nn2;
140 
	if (nn1 <= nn2) {
141 
	    n1 = nn1;
142 
	    n2 = nn2;
143 
	} else {
144 
	    n1 = nn2;
145 
	    n2 = nn1;
146 
	}
147 
    }
148 
    if (setup2) {
149 
	ks = kk;
150 
	if (kk + kk >= tn) {
59097 ripley 151 
	    k = (int)(tn - kk);
574 ihaka 152 
	} else {
153 
	    k = kk;
154 
	}
155 
    }
156 
    if (setup1 || setup2) {
70052 maechler 157 
	m = (int) ((k + 1.) * (n1 + 1.) / (tn + 2.));
574 ihaka 158 
	minjx = imax2(0, k - n2);
159 
	maxjx = imin2(n1, k);
70052 maechler 160 
#ifdef DEBUG_rhyper
161 
	REprintf("rhyper(nn1=%d, nn2=%d, kk=%d), setup: floor(mean)= m=%d, jx in (%d..%d)\n",
162 
		 nn1, nn2, kk, m, minjx, maxjx);
163 
#endif
574 ihaka 164 
    }
4562 pd 165 
    /* generate random variate --- Three basic cases */
574 ihaka 166 
 
4562 pd 167 
    if (minjx == maxjx) { /* I: degenerate distribution ---------------- */
70052 maechler 168 
#ifdef DEBUG_rhyper
169 
	REprintf("rhyper(), branch I (degenerate)\n");
170 
#endif
574 ihaka 171 
	ix = maxjx;
70052 maechler 172 
	goto L_finis; // return appropriate variate
1362 maechler 173 
 
68749 maechler 174 
    } else if (m - minjx < 10) { // II: (Scaled) algorithm HIN (inverse transformation) ----
175 
	const static double scale = 1e25; // scaling factor against (early) underflow
176 
	const static double con = 57.5646273248511421;
68928 ripley 177 
					  // 25*log(10) = log(scale) { <==> exp(con) == scale }
574 ihaka 178 
	if (setup1 || setup2) {
68749 maechler 179 
	    double lw; // log(w);  w = exp(lw) * scale = exp(lw + log(scale)) = exp(lw + con)
574 ihaka 180 
	    if (k < n2) {
68749 maechler 181 
		lw = afc(n2) + afc(n1 + n2 - k) - afc(n2 - k) - afc(n1 + n2);
574 ihaka 182 
	    } else {
68749 maechler 183 
		lw = afc(n1) + afc(     k     ) - afc(k - n2) - afc(n1 + n2);
574 ihaka 184 
	    }
68749 maechler 185 
	    w = exp(lw + con);
574 ihaka 186 
	}
68749 maechler 187 
	double p, u;
70052 maechler 188 
#ifdef DEBUG_rhyper
189 
	REprintf("rhyper(), branch II; w = %g > 0\n", w);
190 
#endif
574 ihaka 191 
      L10:
192 
	p = w;
193 
	ix = minjx;
7725 ripley 194 
	u = unif_rand() * scale;
70052 maechler 195 
#ifdef DEBUG_rhyper
196 
	REprintf("  _new_ u = %g\n", u);
197 
#endif
68749 maechler 198 
	while (u > p) {
13434 maechler 199 
	    u -= p;
68749 maechler 200 
	    p *= ((double) n1 - ix) * (k - ix);
13434 maechler 201 
	    ix++;
574 ihaka 202 
	    p = p / ix / (n2 - k + ix);
70052 maechler 203 
#ifdef DEBUG_rhyper
204 
	    REprintf("       ix=%3d, u=%11g, p=%20.14g (u-p=%g)\n", ix, u, p, u-p);
205 
#endif
574 ihaka 206 
	    if (ix > maxjx)
207 
		goto L10;
70052 maechler 208 
	    // FIXME  if(p == 0.)  we also "have lost"  => goto L10
574 ihaka 209 
	}
68749 maechler 210 
    } else { /* III : H2PE Algorithm --------------------------------------- */
574 ihaka 211 
 
68749 maechler 212 
	double u,v;
213 
 
574 ihaka 214 
	if (setup1 || setup2) {
215 
	    s = sqrt((tn - k) * k * n1 * n2 / (tn - 1) / tn / tn);
216 
 
217 
	    /* remark: d is defined in reference without int. */
218 
	    /* the truncation centers the cell boundaries at 0.5 */
219 
 
220 
	    d = (int) (1.5 * s) + .5;
221 
	    xl = m - d + .5;
222 
	    xr = m + d + .5;
13434 maechler 223 
	    a = afc(m) + afc(n1 - m) + afc(k - m) + afc(n2 - k + m);
574 ihaka 224 
	    kl = exp(a - afc((int) (xl)) - afc((int) (n1 - xl))
225 
		     - afc((int) (k - xl))
226 
		     - afc((int) (n2 - k + xl)));
227 
	    kr = exp(a - afc((int) (xr - 1))
228 
		     - afc((int) (n1 - xr + 1))
229 
		     - afc((int) (k - xr + 1))
230 
		     - afc((int) (n2 - k + xr - 1)));
13434 maechler 231 
	    lamdl = -log(xl * (n2 - k + xl) / (n1 - xl + 1) / (k - xl + 1));
232 
	    lamdr = -log((n1 - xr + 1) * (k - xr + 1) / xr / (n2 - k + xr));
574 ihaka 233 
	    p1 = d + d;
234 
	    p2 = p1 + kl / lamdl;
235 
	    p3 = p2 + kr / lamdr;
236 
	}
70052 maechler 237 
#ifdef DEBUG_rhyper
238 
	REprintf("rhyper(), branch III {accept/reject}: (xl,xr)= (%g,%g); (lamdl,lamdr)= (%g,%g)\n",
239 
		 xl, xr, lamdl,lamdr);
240 
	REprintf("-------- p123= c(%g,%g,%g)\n", p1,p2, p3);
241 
#endif
242 
	int n_uv = 0;
574 ihaka 243 
      L30:
7725 ripley 244 
	u = unif_rand() * p3;
245 
	v = unif_rand();
70052 maechler 246 
	n_uv++;
247 
	if(n_uv >= 10000) {
248 
	    REprintf("rhyper() branch III: giving up after %d rejections", n_uv);
249 
	    ML_ERR_return_NAN;
250 
        }
251 
#ifdef DEBUG_rhyper
252 
	REprintf(" ... L30: new (u=%g, v ~ U[0,1])[%d]\n", u, n_uv);
253 
#endif
254 
 
4562 pd 255 
	if (u < p1) {		/* rectangular region */
59170 ripley 256 
	    ix = (int) (xl + u);
4562 pd 257 
	} else if (u <= p2) {	/* left tail */
59170 ripley 258 
	    ix = (int) (xl + log(v) / lamdl);
574 ihaka 259 
	    if (ix < minjx)
260 
		goto L30;
261 
	    v = v * (u - p1) * lamdl;
4562 pd 262 
	} else {		/* right tail */
59170 ripley 263 
	    ix = (int) (xr - log(v) / lamdr);
574 ihaka 264 
	    if (ix > maxjx)
265 
		goto L30;
266 
	    v = v * (u - p2) * lamdr;
267 
	}
268 
 
269 
	/* acceptance/rejection test */
70052 maechler 270 
	Rboolean reject = TRUE;
574 ihaka 271 
 
272 
	if (m < 100 || ix <= 50) {
273 
	    /* explicit evaluation */
33718 maechler 274 
	    /* The original algorithm (and TOMS 668) have
275 
		   f = f * i * (n2 - k + i) / (n1 - i) / (k - i);
276 
	       in the (m > ix) case, but the definition of the
277 
	       recurrence relation on p134 shows that the +1 is
278 
	       needed. */
70052 maechler 279 
	    int i;
280 
	    double f = 1.0;
574 ihaka 281 
	    if (m < ix) {
282 
		for (i = m + 1; i <= ix; i++)
13434 maechler 283 
		    f = f * (n1 - i + 1) * (k - i + 1) / (n2 - k + i) / i;
574 ihaka 284 
	    } else if (m > ix) {
285 
		for (i = ix + 1; i <= m; i++)
33687 tlumley 286 
		    f = f * i * (n2 - k + i) / (n1 - i + 1) / (k - i + 1);
574 ihaka 287 
	    }
288 
	    if (v <= f) {
10470 maechler 289 
		reject = FALSE;
574 ihaka 290 
	    }
291 
	} else {
70052 maechler 292 
 
293 
	    const static double deltal = 0.0078;
294 
	    const static double deltau = 0.0034;
295 
 
296 
	    double e, g, r, t, y;
68749 maechler 297 
	    double de, dg, dr, ds, dt, gl, gu, nk, nm, ub;
298 
	    double xk, xm, xn, y1, ym, yn, yk, alv;
70052 maechler 299 
 
300 
#ifdef DEBUG_rhyper
301 
	    REprintf(" ... accept/reject 'large' case v=%g\n", v);
302 
#endif
574 ihaka 303 
	    /* squeeze using upper and lower bounds */
304 
	    y = ix;
305 
	    y1 = y + 1.0;
306 
	    ym = y - m;
307 
	    yn = n1 - y + 1.0;
308 
	    yk = k - y + 1.0;
309 
	    nk = n2 - k + y1;
310 
	    r = -ym / y1;
311 
	    s = ym / yn;
312 
	    t = ym / yk;
313 
	    e = -ym / nk;
314 
	    g = yn * yk / (y1 * nk) - 1.0;
315 
	    dg = 1.0;
316 
	    if (g < 0.0)
317 
		dg = 1.0 + g;
318 
	    gu = g * (1.0 + g * (-0.5 + g / 3.0));
319 
	    gl = gu - .25 * (g * g * g * g) / dg;
320 
	    xm = m + 0.5;
321 
	    xn = n1 - m + 0.5;
322 
	    xk = k - m + 0.5;
323 
	    nm = n2 - k + xm;
324 
	    ub = y * gu - m * gl + deltau
325 
		+ xm * r * (1. + r * (-0.5 + r / 3.0))
326 
		+ xn * s * (1. + s * (-0.5 + s / 3.0))
327 
		+ xk * t * (1. + t * (-0.5 + t / 3.0))
328 
		+ nm * e * (1. + e * (-0.5 + e / 3.0));
329 
	    /* test against upper bound */
330 
	    alv = log(v);
331 
	    if (alv > ub) {
10470 maechler 332 
		reject = TRUE;
574 ihaka 333 
	    } else {
334 
				/* test against lower bound */
335 
		dr = xm * (r * r * r * r);
336 
		if (r < 0.0)
13434 maechler 337 
		    dr /= (1.0 + r);
574 ihaka 338 
		ds = xn * (s * s * s * s);
339 
		if (s < 0.0)
13434 maechler 340 
		    ds /= (1.0 + s);
574 ihaka 341 
		dt = xk * (t * t * t * t);
342 
		if (t < 0.0)
13434 maechler 343 
		    dt /= (1.0 + t);
574 ihaka 344 
		de = nm * (e * e * e * e);
345 
		if (e < 0.0)
13434 maechler 346 
		    de /= (1.0 + e);
574 ihaka 347 
		if (alv < ub - 0.25 * (dr + ds + dt + de)
348 
		    + (y + m) * (gl - gu) - deltal) {
10470 maechler 349 
		    reject = FALSE;
33718 maechler 350 
		}
13434 maechler 351 
		else {
352 
		    /* * Stirling's formula to machine accuracy
574 ihaka 353 
		     */
354 
		    if (alv <= (a - afc(ix) - afc(n1 - ix)
355 
				- afc(k - ix) - afc(n2 - k + ix))) {
10470 maechler 356 
			reject = FALSE;
574 ihaka 357 
		    } else {
10470 maechler 358 
			reject = TRUE;
574 ihaka 359 
		    }
360 
		}
361 
	    }
68749 maechler 362 
	} // else
574 ihaka 363 
	if (reject)
364 
	    goto L30;
365 
    }
366 
 
70052 maechler 367 
 
368 
L_finis:
574 ihaka 369 
    /* return appropriate variate */
370 
 
371 
    if (kk + kk >= tn) {
372 
	if (nn1 > nn2) {
373 
	    ix = kk - nn2 + ix;
374 
	} else {
375 
	    ix = nn1 - ix;
376 
	}
377 
    } else {
378 
	if (nn1 > nn2)
379 
	    ix = kk - ix;
380 
    }
381 
    return ix;
382 
}