Skip to content

Commit

Permalink
s/Langage/Language/
Browse files Browse the repository at this point in the history 
git-svn-id: https://svn.r-project.org/R/trunk@4666 00db46b3-68df-0310-9c12-caf00c1e9a41
  • Loading branch information
maechler committed May 21, 1999
1 parent 4ba091e commit 7198da8
Show file tree
Hide file tree
Showing 2 changed files with 115 additions and 111 deletions.
58 changes: 31 additions & 27 deletions src/nmath/qnchisq.c
View file
@@ -1,5 +1,5 @@
/*
* R : A Computer Langage for Statistical Data Analysis
* R : A Computer Language for Statistical Data Analysis
* Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
*
* This program is free software; you can redistribute it and/or modify
Expand All @@ -21,33 +21,37 @@

double qnchisq(double p, double n, double lambda)
{
double ux, lx, nx;
double acu = 1.0e-12;
double ux, lx, nx;
double acu = 1.0e-12;

#ifdef IEEE_754
if (ISNAN(p) || ISNAN(n) || ISNAN(lambda))
return p + n + lambda;
if (!FINITE(n)) {
ML_ERROR(ME_DOMAIN);
return ML_NAN;
}
if (ISNAN(p) || ISNAN(n) || ISNAN(lambda))
return p + n + lambda;
if (!FINITE(n)) {
ML_ERROR(ME_DOMAIN);
return ML_NAN;
}
#endif
n = floor(n + 0.5);
if (p < 0 || p >= 1 || n < 1 || lambda < 0) {
ML_ERROR(ME_DOMAIN);
return ML_NAN;
}
if (p == 0)
return 0;
for (ux = 1.0; pnchisq(ux, n, lambda) < p; ux *= 2);
for (lx = ux; pnchisq(lx, n, lambda) > p; lx *= 0.5);
do {
nx = 0.5 * (lx + ux);
if (pnchisq(nx, n, lambda) > p)
ux = nx;
else
lx = nx;
}
while ((ux - lx) / nx > acu);
return 0.5 * (ux + lx);
n = floor(n + 0.5);
if (p < 0 || p >= 1 || n < 1 || lambda < 0) {
ML_ERROR(ME_DOMAIN);
return ML_NAN;
}
if (p == 0)
return 0;
/* Invert pnchisq(.) finding an upper and lower bound;
then interval halfing : */

/* FIXME: Use less precision here (see pnchisq() !) */
for (ux = 1.0; pnchisq(ux, n, lambda) < p; ux *= 2);
for (lx = ux; pnchisq(lx, n, lambda) > p; lx *= 0.5);
do {
nx = 0.5 * (lx + ux);
if (pnchisq(nx, n, lambda) > p)
ux = nx;
else
lx = nx;
}
while ((ux - lx) / nx > acu);
return 0.5 * (ux + lx);
}
168 changes: 84 additions & 84 deletions src/nmath/rbeta.c
View file
@@ -1,5 +1,5 @@
/*
* R : A Computer Langage for Statistical Data Analysis
* R : A Computer Language for Statistical Data Analysis
* Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
*
* This program is free software; you can redistribute it and/or modify
Expand Down Expand Up @@ -32,94 +32,94 @@ static double expmax = 0.0;

double rbeta(double aa, double bb)
{
static double a, b, delta, r, s, t, u1, u2, v, w, y, z;
static double alpha, beta, gamma, k1, k2;
static double olda = -1.0;
static double oldb = -1.0;
int qsame;
static double a, b, delta, r, s, t, u1, u2, v, w, y, z;
static double alpha, beta, gamma, k1, k2;
static double olda = -1.0;
static double oldb = -1.0;
int qsame;

if (expmax == 0.0)
expmax = log(DBL_MAX);
if (expmax == 0.0)
expmax = log(DBL_MAX);

qsame = (olda == aa) && (oldb == bb);
qsame = (olda == aa) && (oldb == bb);

if (!qsame) {
if (aa > 0.0 && bb > 0.0) {
olda = aa;
oldb = bb;
} else {
ML_ERROR(ME_DOMAIN);
return ML_NAN;
}
}
if (fmin2(aa, bb) <= 1.0) { /* Algorithm BC */
if (!qsame) {
a = fmax2(aa, bb);
b = fmin2(aa, bb);
alpha = a + b;
beta = 1.0 / b;
delta = 1.0 + a - b;
k1 = delta * (0.0138889 + 0.0416667 * b) /
(a * beta - 0.777778);
k2 = 0.25 + (0.5 + 0.25 / delta) * b;
}
repeat {
u1 = sunif();
u2 = sunif();
if (u1 < 0.5) {
y = u1 * u2;
z = u1 * y;
if (0.25 * u2 + z - y >= k1)
continue;
} else {
z = u1 * u1 * u2;
if (z <= 0.25)
break;
if (z >= k2)
continue;
}
v = beta * log(u1 / (1.0 - u1));
if (v <= expmax)
w = a * exp(v);
else
w = DBL_MAX;
if (alpha * (log(alpha / (b + w)) + v) - 1.3862944
>= log(z))
goto deliver;
}
v = beta * log(u1 / (1.0 - u1));
if (v <= expmax)
w = a * exp(v);
else
w = DBL_MAX;
} else { /* Algorithm BB */
if (!qsame) {
if (aa > 0.0 && bb > 0.0) {
olda = aa;
oldb = bb;
} else {
ML_ERROR(ME_DOMAIN);
return ML_NAN;
}
a = fmin2(aa, bb);
b = fmax2(aa, bb);
alpha = a + b;
beta = sqrt((alpha - 2.0) / (2.0 * a * b - alpha));
gamma = a + 1.0 / beta;
}
if (fmin2(aa, bb) <= 1.0) { /* Algorithm BC */
if (!qsame) {
a = fmax2(aa, bb);
b = fmin2(aa, bb);
alpha = a + b;
beta = 1.0 / b;
delta = 1.0 + a - b;
k1 = delta * (0.0138889 + 0.0416667 * b) /
(a * beta - 0.777778);
k2 = 0.25 + (0.5 + 0.25 / delta) * b;
}
repeat {
u1 = sunif();
u2 = sunif();
if (u1 < 0.5) {
y = u1 * u2;
z = u1 * y;
if (0.25 * u2 + z - y >= k1)
continue;
} else {
z = u1 * u1 * u2;
if (z <= 0.25)
break;
if (z >= k2)
continue;
}
v = beta * log(u1 / (1.0 - u1));
if (v <= expmax)
w = a * exp(v);
else
w = DBL_MAX;
if (alpha * (log(alpha / (b + w)) + v) - 1.3862944
>= log(z))
goto deliver;
}
v = beta * log(u1 / (1.0 - u1));
if (v <= expmax)
w = a * exp(v);
else
w = DBL_MAX;
} else { /* Algorithm BB */
if (!qsame) {
a = fmin2(aa, bb);
b = fmax2(aa, bb);
alpha = a + b;
beta = sqrt((alpha - 2.0) / (2.0 * a * b - alpha));
gamma = a + 1.0 / beta;
}
do {
u1 = sunif();
u2 = sunif();
v = beta * log(u1 / (1.0 - u1));
if (v <= expmax)
w = a * exp(v);
else
w = DBL_MAX;
z = u1 * u1 * u2;
r = gamma * v - 1.3862944;
s = a + r - w;
if (s + 2.609438 >= 5.0 * z)
break;
t = log(z);
if (s > t)
break;
}
while (r + alpha * log(alpha / (b + w)) < t);
do {
u1 = sunif();
u2 = sunif();
v = beta * log(u1 / (1.0 - u1));
if (v <= expmax)
w = a * exp(v);
else
w = DBL_MAX;
z = u1 * u1 * u2;
r = gamma * v - 1.3862944;
s = a + r - w;
if (s + 2.609438 >= 5.0 * z)
break;
t = log(z);
if (s > t)
break;
}
while (r + alpha * log(alpha / (b + w)) < t);
}

deliver:
return (aa != a) ? b / (b + w) : w / (b + w);
deliver:
return (aa != a) ? b / (b + w) : w / (b + w);
}

0 comments on commit 7198da8

Please sign in to comment.