The R Project SVN R

Rev

Rev 59039 | Rev 64668 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 59039 Rev 59097
Line 52... Line 52...
52 
    static int nsave = -1;
 52 
    static int nsave = -1;
53 
    static int m;
 53 
    static int m;
54 
 
 54 
 
55 
    double f, f1, f2, u, v, w, w2, x, x1, x2, z, z2;
 55 
    double f, f1, f2, u, v, w, w2, x, x1, x2, z, z2;
56 
    double p, q, np, g, r, al, alv, amaxp, ffm, ynorm;
 56 
    double p, q, np, g, r, al, alv, amaxp, ffm, ynorm;
57 
    int i,ix,k, n;
 57 
    int i, ix, k, n;
58 
 
 58 
 
59 
    if (!R_FINITE(nin)) ML_ERR_return_NAN;
 59 
    if (!R_FINITE(nin)) ML_ERR_return_NAN;
60 
    r = floor(nin + 0.5);
 60 
    r = floor(nin + 0.5);
61 
    if (r != nin) ML_ERR_return_NAN;
 61 
    if (r != nin) ML_ERR_return_NAN;
62 
    if (!R_FINITE(pp) ||
 62 
    if (!R_FINITE(pp) ||
Line 68... Line 68...
68 
 
 68 
 
69 
    if (r >= INT_MAX)/* evade integer overflow,
 69 
    if (r >= INT_MAX)/* evade integer overflow,
70 
			and r == INT_MAX gave only even values */
 70 
			and r == INT_MAX gave only even values */
71 
	return qbinom(unif_rand(), r, pp, /*lower_tail*/ 0, /*log_p*/ 0);
 71 
	return qbinom(unif_rand(), r, pp, /*lower_tail*/ 0, /*log_p*/ 0);
72 
    /* else */
 72 
    /* else */
73 
    n = r;
 73 
    n = (int) r;
74 
 
 74 
 
75 
    p = fmin2(pp, 1. - pp);
 75 
    p = fmin2(pp, 1. - pp);
76 
    q = 1. - p;
 76 
    q = 1. - p;
77 
    np = n * p;
 77 
    np = n * p;
78 
    r = p / q;
 78 
    r = p / q;
Line 90... Line 90...
90 
	    /* inverse cdf logic for mean less than 30 */
 90 
	    /* inverse cdf logic for mean less than 30 */
91 
	    qn = pow(q, (double) n);
 91 
	    qn = pow(q, (double) n);
92 
	    goto L_np_small;
 92 
	    goto L_np_small;
93 
	} else {
 93 
	} else {
94 
	    ffm = np + p;
 94 
	    ffm = np + p;
95 
	    m = ffm;
 95 
	    m = (int) ffm;
96 
	    fm = m;
 96 
	    fm = m;
97 
	    npq = np * q;
 97 
	    npq = np * q;
98 
	    p1 = (int)(2.195 * sqrt(npq) - 4.6 * q) + 0.5;
 98 
	    p1 = (int)(2.195 * sqrt(npq) - 4.6 * q) + 0.5;
99 
	    xm = fm + 0.5;
 99 
	    xm = fm + 0.5;
100 
	    xl = xm - p1;
 100 
	    xl = xm - p1;
Line 117... Line 117...
117 
    repeat {
 117 
    repeat {
118 
      u = unif_rand() * p4;
 118 
      u = unif_rand() * p4;
119 
      v = unif_rand();
 119 
      v = unif_rand();
120 
      /* triangular region */
 120 
      /* triangular region */
121 
      if (u <= p1) {
 121 
      if (u <= p1) {
122 
	  ix = xm - p1 * v + u;
 122 
	  ix = (int)(xm - p1 * v + u);
123 
	  goto finis;
 123 
	  goto finis;
124 
      }
 124 
      }
125 
      /* parallelogram region */
 125 
      /* parallelogram region */
126 
      if (u <= p2) {
 126 
      if (u <= p2) {
127 
	  x = xl + (u - p1) / c;
 127 
	  x = xl + (u - p1) / c;
128 
	  v = v * c + 1.0 - fabs(xm - x) / p1;
 128 
	  v = v * c + 1.0 - fabs(xm - x) / p1;
129 
	  if (v > 1.0 || v <= 0.)
 129 
	  if (v > 1.0 || v <= 0.)
130 
	      continue;
 130 
	      continue;
131 
	  ix = x;
 131 
	  ix = (int) x;
132 
      } else {
 132 
      } else {
133 
	  if (u > p3) {	/* right tail */
 133 
	  if (u > p3) {	/* right tail */
134 
	      ix = xr - log(v) / xlr;
 134 
	      ix = (int)(xr - log(v) / xlr);
135 
	      if (ix > n)
 135 
	      if (ix > n)
136 
		  continue;
 136 
		  continue;
137 
	      v = v * (u - p3) * xlr;
 137 
	      v = v * (u - p3) * xlr;
138 
	  } else {/* left tail */
 138 
	  } else {/* left tail */
139 
	      ix = xl + log(v) / xll;
 139 
	      ix = (int)(xl + log(v) / xll);
140 
	      if (ix < 0)
 140 
	      if (ix < 0)
141 
		  continue;
 141 
		  continue;
142 
	      v = v * (u - p2) * xll;
 142 
	      v = v * (u - p2) * xll;
143 
	  }
 143 
	  }
144 
      }
 144 
      }