Commit 256affb8 authored by Morten Welinder's avatar Morten Welinder
Browse files

Massive update for statistical functions.

Using code from "R" gives better precision and we avoid a lot of root-
finding loops, i.e., we get speed too.
parent 3e62f1f3
1999-05-25 Morten Welinder <terra@diku.dk>
* src/mathfunc.c, src/mathfunc.h: New files, mostly taken from the
R package. (It's a GPL'ed gold mine.)
* src/*.c: Use <math.h>, not "math.h".
* src/utils.c (random_normal): Use qnorm, not inv_phi.
* src/fn-stat.c: Move all R code to mathfunc.c
Change all uses to phi to pnorm.
(gnumeric_normsinv): Use qnorm.
(gnumeric_confidence): Use qnorm.
(normsinv): Superseded by qnorm.
(help_lognormdist): Fix.
(gnumeric_lognormdist): Fix domain.
(gnumeric_lognormdist): Use plnorm.
(gnumeric_loginv): Use qlnorm.
(gnumeric_norminv): Use qnorm.
(gnumeric_tinv): Use qt.
(gnumeric_fdist): Use qf.
(gnumeric_gammainv): Use qgamma.
(gnumeric_chiinv): Use qchisq.
* src/Makefile.am (GNUMERIC_BASE_SOURCES): Add mathfunc.c and
mathfunc.h.
1999-05-25 Morten Welinder <terra@diku.dk>
* src/fn-string.c (gnumeric_code): Handle compilers for which the
......
1999-05-25 Morten Welinder <terra@diku.dk>
* src/mathfunc.c, src/mathfunc.h: New files, mostly taken from the
R package. (It's a GPL'ed gold mine.)
* src/*.c: Use <math.h>, not "math.h".
* src/utils.c (random_normal): Use qnorm, not inv_phi.
* src/fn-stat.c: Move all R code to mathfunc.c
Change all uses to phi to pnorm.
(gnumeric_normsinv): Use qnorm.
(gnumeric_confidence): Use qnorm.
(normsinv): Superseded by qnorm.
(help_lognormdist): Fix.
(gnumeric_lognormdist): Fix domain.
(gnumeric_lognormdist): Use plnorm.
(gnumeric_loginv): Use qlnorm.
(gnumeric_norminv): Use qnorm.
(gnumeric_tinv): Use qt.
(gnumeric_fdist): Use qf.
(gnumeric_gammainv): Use qgamma.
(gnumeric_chiinv): Use qchisq.
* src/Makefile.am (GNUMERIC_BASE_SOURCES): Add mathfunc.c and
mathfunc.h.
1999-05-25 Morten Welinder <terra@diku.dk>
* src/fn-string.c (gnumeric_code): Handle compilers for which the
......
1999-05-25 Morten Welinder <terra@diku.dk>
* src/mathfunc.c, src/mathfunc.h: New files, mostly taken from the
R package. (It's a GPL'ed gold mine.)
* src/*.c: Use <math.h>, not "math.h".
* src/utils.c (random_normal): Use qnorm, not inv_phi.
* src/fn-stat.c: Move all R code to mathfunc.c
Change all uses to phi to pnorm.
(gnumeric_normsinv): Use qnorm.
(gnumeric_confidence): Use qnorm.
(normsinv): Superseded by qnorm.
(help_lognormdist): Fix.
(gnumeric_lognormdist): Fix domain.
(gnumeric_lognormdist): Use plnorm.
(gnumeric_loginv): Use qlnorm.
(gnumeric_norminv): Use qnorm.
(gnumeric_tinv): Use qt.
(gnumeric_fdist): Use qf.
(gnumeric_gammainv): Use qgamma.
(gnumeric_chiinv): Use qchisq.
* src/Makefile.am (GNUMERIC_BASE_SOURCES): Add mathfunc.c and
mathfunc.h.
1999-05-25 Morten Welinder <terra@diku.dk>
* src/fn-string.c (gnumeric_code): Handle compilers for which the
......
1999-05-25 Morten Welinder <terra@diku.dk>
* src/mathfunc.c, src/mathfunc.h: New files, mostly taken from the
R package. (It's a GPL'ed gold mine.)
* src/*.c: Use <math.h>, not "math.h".
* src/utils.c (random_normal): Use qnorm, not inv_phi.
* src/fn-stat.c: Move all R code to mathfunc.c
Change all uses to phi to pnorm.
(gnumeric_normsinv): Use qnorm.
(gnumeric_confidence): Use qnorm.
(normsinv): Superseded by qnorm.
(help_lognormdist): Fix.
(gnumeric_lognormdist): Fix domain.
(gnumeric_lognormdist): Use plnorm.
(gnumeric_loginv): Use qlnorm.
(gnumeric_norminv): Use qnorm.
(gnumeric_tinv): Use qt.
(gnumeric_fdist): Use qf.
(gnumeric_gammainv): Use qgamma.
(gnumeric_chiinv): Use qchisq.
* src/Makefile.am (GNUMERIC_BASE_SOURCES): Add mathfunc.c and
mathfunc.h.
1999-05-25 Morten Welinder <terra@diku.dk>
* src/fn-string.c (gnumeric_code): Handle compilers for which the
......
......@@ -6,7 +6,7 @@
*/
#include <config.h>
#include <gnome.h>
#include "math.h"
#include <math.h>
#include "gnumeric.h"
#include "gnumeric-sheet.h"
#include "utils.h"
......
......@@ -7,7 +7,7 @@
*/
#include <config.h>
#include <gnome.h>
#include "math.h"
#include <math.h>
#include "gnumeric.h"
#include "gnumeric-sheet.h"
#include "utils.h"
......
......@@ -8,7 +8,7 @@
#include <config.h>
#include <gnome.h>
#include <ctype.h>
#include "math.h"
#include <math.h>
#include "numbers.h"
#include "gnumeric.h"
#include "gnumeric-sheet.h"
......
......@@ -9,7 +9,7 @@
*/
#include <config.h>
#include <gnome.h>
#include "math.h"
#include <math.h>
#include "gnumeric.h"
#include "gnumeric-sheet.h"
#include "utils.h"
......
......@@ -8,7 +8,7 @@
#include <config.h>
#include <gnome.h>
#include <ctype.h>
#include "math.h"
#include <math.h>
#include "gnumeric.h"
#include "gnumeric-sheet.h"
#include "utils.h"
......
......@@ -7,7 +7,6 @@
*/
#include <config.h>
#include <gnome.h>
#include "math.h"
#include "gnumeric.h"
#include "gnumeric-sheet.h"
#include "utils.h"
......
......@@ -7,7 +7,7 @@
#include <config.h>
#include <gnome.h>
#include "math.h"
#include <math.h>
#include "numbers.h"
#include "gnumeric.h"
#include "gnumeric-sheet.h"
......
......@@ -8,23 +8,14 @@
*/
#include <config.h>
#include <gnome.h>
#include "math.h"
#include <math.h>
#include "mathfunc.h"
#include "numbers.h"
#include "gnumeric.h"
#include "gnumeric-sheet.h"
#include "utils.h"
#include "func.h"
#define M_LN_SQRT_2PI 0.918938533204672741780329736406 /* log(sqrt(2*pi)) */
/* The cumulative distribution function for the 0-1 normal distribution. */
static float_t
phi (float_t x)
{
return (1.0 + erf (x / M_SQRT2)) / 2.0;
}
static guint
float_hash (const float_t *d)
{
......@@ -61,448 +52,6 @@ float_compare_d (const float_t *a, const float_t *b)
return 1;
}
#define fmin2(a,b) MIN(a,b)
#define fmax2(a,b) MAX(a,b)
/* This function is originally taken from R package
* (src/nmath/pgamma.c) written and copyrighted (1998) by Ross Ihaka.
* Modified for Gnumeric by Jukka-Pekka Iivonen
*/
static float_t
pgamma(double x, double p, double scale)
{
const float_t third = 1.0 / 3.0;
const float_t xbig = 1.0e+8;
const float_t oflo = 1.0e+37;
const float_t plimit = 1000.0e0;
const float_t elimit = -88.0e0;
float_t pn1, pn2, pn3, pn4, pn5, pn6, arg, c, rn, a, b, an;
float_t sum;
x = x / scale;
if (x <= 0)
return 0.0;
/* use a normal approximation if p > plimit */
if (p > plimit) {
pn1 = sqrt(p) * 3.0 * (pow(x/p, third) + 1.0 /
(p * 9.0) - 1.0);
return phi(pn1);
}
/* if x is extremely large compared to p then return 1 */
if (x > xbig)
return 1.0;
if (x <= 1.0 || x < p) {
/* use pearson's series expansion. */
arg = p * log(x) - x - lgamma(p + 1.0);
c = 1.0;
sum = 1.0;
a = p;
do {
a = a + 1.0;
c = c * x / a;
sum = sum + c;
} while (c > DBL_EPSILON);
arg = arg + log(sum);
sum = 0;
if (arg >= elimit)
sum = exp(arg);
} else {
/* use a continued fraction expansion */
arg = p * log(x) - x - lgamma(p);
a = 1.0 - p;
b = a + x + 1.0;
c = 0;
pn1 = 1.0;
pn2 = x;
pn3 = x + 1.0;
pn4 = x * b;
sum = pn3 / pn4;
for (;;) {
a = a + 1.0;
b = b + 2.0;
c = c + 1.0;
an = a * c;
pn5 = b * pn3 - an * pn1;
pn6 = b * pn4 - an * pn2;
if (fabs(pn6) > 0) {
rn = pn5 / pn6;
if (fabs(sum - rn) <=
fmin2(DBL_EPSILON, DBL_EPSILON * rn))
break;
sum = rn;
}
pn1 = pn3;
pn2 = pn4;
pn3 = pn5;
pn4 = pn6;
if (fabs(pn5) >= oflo) {
/* re-scale the terms in continued fraction */
/* if they are large */
pn1 = pn1 / oflo;
pn2 = pn2 / oflo;
pn3 = pn3 / oflo;
pn4 = pn4 / oflo;
}
}
arg = arg + log(sum);
sum = 1.0;
if (arg >= elimit)
sum = 1.0 - exp(arg);
}
return sum;
}
/* This function is originally taken from R package
* (src/nmath/i1match.c) written and copyrighted (1998) by Ross Ihaka.
* Modified for Gnumeric by Jukka-Pekka Iivonen
*/
static int
i1mach (int i)
{
switch(i) {
case 1:
return 5;
case 2:
return 6;
case 3:
return 0;
case 4:
return 0;
case 5:
return CHAR_BIT * sizeof(int);
case 6:
return sizeof(int)/sizeof(char);
case 7:
return 2;
case 8:
return CHAR_BIT * sizeof(int) - 1;
case 9:
return INT_MAX;
case 10:
return FLT_RADIX;
case 11:
return FLT_MANT_DIG;
case 12:
return FLT_MIN_EXP;
case 13:
return FLT_MAX_EXP;
case 14:
return DBL_MANT_DIG;
case 15:
return DBL_MIN_EXP;
case 16:
return DBL_MAX_EXP;
default:
return 0;
}
}
/* This function is originally taken from R package
* (src/nmath/d1match.c) written and copyrighted (1998) by Ross Ihaka.
* Modified for Gnumeric by Jukka-Pekka Iivonen
*/
static float_t
d1mach(int i)
{
switch(i) {
case 1:
return DBL_MIN;
case 2:
return DBL_MAX;
case 3:
return pow((double)i1mach(10), -(double)i1mach(14));
case 4:
return pow((double)i1mach(10), 1-(double)i1mach(14));
case 5:
return log10(2.0);
default:
return 0.0;
}
}
/* This function is originally taken from R package
* (src/nmath/chebyshev.c) written and copyrighted (1998) by Ross Ihaka.
* Modified for Gnumeric by Jukka-Pekka Iivonen
*/
static int
chebyshev_init(float_t *dos, int nos, float_t eta)
{
int i, ii;
float_t err;
if (nos < 1)
return 0;
err = 0.0;
i = 0; /* just to avoid compiler warnings */
for (ii=1; ii<=nos; ii++) {
i = nos - ii;
err += fabs(dos[i]);
if (err > eta) {
return i;
}
}
return i;
}
/* This function is originally taken from R package
* (src/nmath/chebyshev.c) written and copyrighted (1998) by Ross Ihaka.
* Modified for Gnumeric by Jukka-Pekka Iivonen
*/
static float_t
chebyshev_eval(float_t x, float_t *a, int n)
{
float_t b0, b1, b2, twox;
int i;
twox = x * 2;
b2 = b1 = 0;
b0 = 0;
for (i = 1; i <= n; i++) {
b2 = b1;
b1 = b0;
b0 = twox * b1 - b2 + a[n - i];
}
return (b0 - b2) * 0.5;
}
/* This function is originally taken from R package
* (src/nmath/lgammacor.c) written and copyrighted (1998) by Ross Ihaka.
* Modified for Gnumeric by Jukka-Pekka Iivonen
*/
static float_t
lgammacor(float_t x)
{
static double algmcs[15] = {
+.1666389480451863247205729650822e+0,
-.1384948176067563840732986059135e-4,
+.9810825646924729426157171547487e-8,
-.1809129475572494194263306266719e-10,
+.6221098041892605227126015543416e-13,
-.3399615005417721944303330599666e-15,
+.2683181998482698748957538846666e-17,
-.2868042435334643284144622399999e-19,
+.3962837061046434803679306666666e-21,
-.6831888753985766870111999999999e-23,
+.1429227355942498147573333333333e-24,
-.3547598158101070547199999999999e-26,
+.1025680058010470912000000000000e-27,
-.3401102254316748799999999999999e-29,
+.1276642195630062933333333333333e-30
};
static int nalgm = 0;
static float_t xbig = 0;
static float_t xmax = 0;
float_t tmp;
if (nalgm == 0) {
nalgm = chebyshev_init(algmcs, 15, d1mach(3));
xbig = 1 / sqrt(d1mach(3));
xmax = exp(fmin2(log(d1mach(2) / 12), -log(12 * d1mach(1))));
}
if (x < xbig) {
tmp = 10 / x;
return chebyshev_eval(tmp * tmp * 2 - 1, algmcs, nalgm) / x;
}
else
return (1 / (x * 12));
}
/* This function is originally taken from R package
* (src/nmath/lbeta.c) written and copyrighted (1998) by Ross Ihaka.
* Modified for Gnumeric by Jukka-Pekka Iivonen
*/
static float_t
lbeta(float_t a, float_t b)
{
static float_t corr, p, q;
p = q = a;
if (b < p)
p = b; /* := min(a,b) */
if (b > q)
q = b; /* := max(a,b) */
/* both arguments must be >= 0 */
if (p >= 10) {
/* p and q are big. */
corr = lgammacor(p) + lgammacor(q) - lgammacor(p + q);
return log(q) * -0.5 + M_LN_SQRT_2PI + corr
+ (p - 0.5) * log(p / (p + q)) + q * log(1+(-p / (p + q)));
}
else if (q >= 10) {
/* p is small, but q is big. */
corr = lgammacor(q) - lgammacor(p + q);
return lgamma(p) + corr + p - p * log(p + q)
+ (q - 0.5) * log(1+(-p / (p + q)));
}
else
/* p and q are small: p <= q > 10. */
return log(exp(lgamma(p)) *
(exp(lgamma(q)) / exp(lgamma(p + q))));
}
/* This function is originally taken from R package
* (src/nmath/pbeta.c) written and copyrighted (1998) by Ross Ihaka.
* Modified for Gnumeric by Jukka-Pekka Iivonen
*/
static float_t
pbeta_raw(float_t x, float_t pin, float_t qin)
{
float_t ans, c, finsum, p, ps, p1, q, term, xb, xi, y;
int n, i, ib;
static float_t eps = 0;
static float_t alneps = 0;
static float_t sml = 0;
static float_t alnsml = 0;
if (eps == 0) { /* initialize machine constants ONCE */
eps = d1mach(3);
alneps = log(eps);
sml = d1mach(1);
alnsml = log(sml);
}
y = x;
p = pin;
q = qin;
/* swap tails if x is greater than the mean */
if (p / (p + q) < x) {
y = 1 - y;
p = qin;
q = pin;
}
if ((p + q) * y / (p + 1) < eps) {
/* tail approximation */
ans = 0;
xb = p * log(fmax2(y, sml)) - log(p) - lbeta(p, q);
if (xb > alnsml && y != 0)
ans = exp(xb);
if (y != x || p != pin)
ans = 1 - ans;
}
else {
/*___ FIXME ___: This takes forever (or ends wrongly)
* when (one or) both p & q are huge
*/
/* evaluate the infinite sum first. term will equal */
/* y^p / beta(ps, p) * (1 - ps)-sub-i * y^i / fac(i) */
ps = q - floor(q);
if (ps == 0)
ps = 1;
xb = p * log(y) - lbeta(ps, p) - log(p);
ans = 0;
if (xb >= alnsml) {
ans = exp(xb);
term = ans * p;
if (ps != 1) {
n = fmax2(alneps/log(y), 4.0);
for(i=1 ; i<= n ; i++) {
xi = i;
term = term * (xi - ps) * y / xi;
ans = ans + term / (p + xi);
}
}
}
/* now evaluate the finite sum, maybe. */
if (q > 1) {
xb = p * log(y) + q * log(1 - y) - lbeta(p, q) - log(q);
ib = fmax2(xb / alnsml, 0.0);
term = exp(xb - ib * alnsml);
c = 1 / (1 - y);
p1 = q * c / (p + q - 1);
finsum = 0;
n = q;
if (q == n)
n = n - 1;
for(i=1 ; i<=n ; i++) {
if (p1 <= 1 && term / eps <= finsum)
break;
xi = i;
term = (q - xi + 1) * c * term / (p + q - xi);
if (term > 1) {
ib = ib - 1;
term = term * sml;
}
if (ib == 0)
finsum = finsum + term;
}
ans = ans + finsum;
}
if (y != x || p != pin)
ans = 1 - ans;
ans = fmax2(fmin2(ans, 1.0), 0.0);
}
return ans;