The R Project SVN R

Rev

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

Rev 3076 Rev 4562
Line 114... Line 114...
114 
    }
 114 
    }
115 
#endif
 115 
#endif
116 
 
 116 
 
117 
    nn1 = floor(nn1in+0.5);
 117 
    nn1 = floor(nn1in+0.5);
118 
    nn2 = floor(nn2in+0.5);
 118 
    nn2 = floor(nn2in+0.5);
119 
    kk = floor(kkin+0.5);
 119 
    kk	= floor(kkin +0.5);
120 
 
 120 
 
121 
    if (nn1 < 0 || nn2 < 0 || kk < 0 || kk > nn1 + nn2) {
 121 
    if (nn1 < 0 || nn2 < 0 || kk < 0 || kk > nn1 + nn2) {
122 
	ML_ERROR(ME_DOMAIN);
 122 
	ML_ERROR(ME_DOMAIN);
123 
	return ML_NAN;
 123 
	return ML_NAN;
124 
    }
 124 
    }
Line 156... Line 156...
156 
    if (setup1 || setup2) {
 156 
    if (setup1 || setup2) {
157 
	m = (k + 1.0) * (n1 + 1.0) / (tn + 2.0);
 157 
	m = (k + 1.0) * (n1 + 1.0) / (tn + 2.0);
158 
	minjx = imax2(0, k - n2);
 158 
	minjx = imax2(0, k - n2);
159 
	maxjx = imin2(n1, k);
 159 
	maxjx = imin2(n1, k);
160 
    }
 160 
    }
161 
    /* generate random variate */
 161 
    /* generate random variate --- Three basic cases */
162 
 
 162 
 
163 
    if (minjx == maxjx) {
 - 
 
164 
	/* degenerate distribution */
 163 
    if (minjx == maxjx) { /* I: degenerate distribution ---------------- */
165 
	ix = maxjx;
 164 
	ix = maxjx;
166 
	/* return ix;
 165 
	/* return ix;
167 
	   No, need to unmangle <TSL>*/
 166 
	   No, need to unmangle <TSL>*/
168 
	/* return appropriate variate */
 167 
	/* return appropriate variate */
169 
 
 168 
 
Line 177... Line 176...
177 
	  if (nn1 > nn2)
 176 
	  if (nn1 > nn2)
178 
	    ix = kk - ix;
 177 
	    ix = kk - ix;
179 
	}
 178 
	}
180 
	return ix;
 179 
	return ix;
181 
 
 180 
 
182 
    } else if (m - minjx < 10) {
 181 
    } else if (m - minjx < 10) { /* II: inverse transformation ---------- */
183 
	/* inverse transformation */
 - 
 
184 
	if (setup1 || setup2) {
 182 
	if (setup1 || setup2) {
185 
	    if (k < n2) {
 183 
	    if (k < n2) {
186 
		w = exp(con + afc(n2) + afc(n1 + n2 - k)
 184 
		w = exp(con + afc(n2) + afc(n1 + n2 - k)
187 
			- afc(n2 - k) - afc(n1 + n2));
 185 
			- afc(n2 - k) - afc(n1 + n2));
188 
	    } else {
 186 
	    } else {
Line 202... Line 200...
202 
	    p = p / ix / (n2 - k + ix);
 200 
	    p = p / ix / (n2 - k + ix);
203 
	    if (ix > maxjx)
 201 
	    if (ix > maxjx)
204 
		goto L10;
 202 
		goto L10;
205 
	    goto L20;
 203 
	    goto L20;
206 
	}
 204 
	}
207 
    } else {
 205 
    } else { /* III : h2pe --------------------------------------------- */
208 
	/* h2pe */
 - 
 
209 
 
 206 
 
210 
	if (setup1 || setup2) {
 207 
	if (setup1 || setup2) {
211 
	    s = sqrt((tn - k) * k * n1 * n2 / (tn - 1) / tn / tn);
 208 
	    s = sqrt((tn - k) * k * n1 * n2 / (tn - 1) / tn / tn);
212 
 
 209 
 
213 
	    /* remark: d is defined in reference without int. */
 210 
	    /* remark: d is defined in reference without int. */
Line 234... Line 231...
234 
	    p3 = p2 + kr / lamdr;
 231 
	    p3 = p2 + kr / lamdr;
235 
	}
 232 
	}
236 
      L30:
 233 
      L30:
237 
	u = sunif() * p3;
 234 
	u = sunif() * p3;
238 
	v = sunif();
 235 
	v = sunif();
239 
	if (u < p1) {
 - 
 
240 
	    /* rectangular region */
 236 
	if (u < p1) {		/* rectangular region */
241 
	    ix = xl + u;
 237 
	    ix = xl + u;
242 
	} else if (u <= p2) {
 238 
	} else if (u <= p2) {	/* left tail */
243 
	    /* left tail */
 - 
 
244 
	    ix = xl + log(v) / lamdl;
 239 
	    ix = xl + log(v) / lamdl;
245 
	    if (ix < minjx)
 240 
	    if (ix < minjx)
246 
		goto L30;
 241 
		goto L30;
247 
	    v = v * (u - p1) * lamdl;
 242 
	    v = v * (u - p1) * lamdl;
248 
	} else {
 - 
 
249 
	    /* right tail */
 243 
	} else {		/* right tail */
250 
	    ix = xr - log(v) / lamdr;
 244 
	    ix = xr - log(v) / lamdr;
251 
	    if (ix > maxjx)
 245 
	    if (ix > maxjx)
252 
		goto L30;
 246 
		goto L30;
253 
	    v = v * (u - p2) * lamdr;
 247 
	    v = v * (u - p2) * lamdr;
254 
	}
 248 
	}