350 lines
8.2 KiB
C
350 lines
8.2 KiB
C
#include "S.h"
|
|
#include "loess.h"
|
|
|
|
/* If your compiler is so ancient it doesn't recognize void, say
|
|
#define void
|
|
*/
|
|
|
|
void
|
|
anova(one, two, out)
|
|
struct loess_struct *one, *two;
|
|
struct anova_struct *out;
|
|
{
|
|
double one_d1, one_d2, one_s, two_d1, two_d2, two_s,
|
|
rssdiff, d1diff, tmp, pf();
|
|
int max_enp;
|
|
|
|
one_d1 = one->out.one_delta;
|
|
one_d2 = one->out.two_delta;
|
|
one_s = one->out.s;
|
|
two_d1 = two->out.one_delta;
|
|
two_d2 = two->out.two_delta;
|
|
two_s = two->out.s;
|
|
|
|
rssdiff = fabs(one_s * one_s * one_d1 - two_s * two_s * two_d1);
|
|
d1diff = fabs(one_d1 - two_d1);
|
|
out->dfn = d1diff * d1diff / fabs(one_d2 - two_d2);
|
|
max_enp = (one->out.enp > two->out.enp);
|
|
tmp = max_enp ? one_d1 : two_d1;
|
|
out->dfd = tmp * tmp / (max_enp ? one_d2 : two_d2);
|
|
tmp = max_enp ? one_s : two_s;
|
|
out->F_value = (rssdiff / d1diff) / (tmp * tmp);
|
|
out->Pr_F = 1 - pf(out->F_value, out->dfn, out->dfd);
|
|
}
|
|
|
|
void
|
|
pointwise(pre, m, coverage, ci)
|
|
struct pred_struct *pre;
|
|
int m;
|
|
double coverage;
|
|
struct ci_struct *ci;
|
|
{
|
|
double t_dist, limit, fit, qt();
|
|
int i;
|
|
|
|
ci->fit = (double *) malloc(m * sizeof(double));
|
|
ci->upper = (double *) malloc(m * sizeof(double));
|
|
ci->lower = (double *) malloc(m * sizeof(double));
|
|
|
|
t_dist = qt(1 - (1 - coverage)/2, pre->df);
|
|
for(i = 0; i < m; i++) {
|
|
limit = pre->se_fit[i] * t_dist;
|
|
ci->fit[i] = fit = pre->fit[i];
|
|
ci->upper[i] = fit + limit;
|
|
ci->lower[i] = fit - limit;
|
|
}
|
|
}
|
|
|
|
void
|
|
pw_free_mem(ci)
|
|
struct ci_struct *ci;
|
|
{
|
|
free(ci->fit);
|
|
free(ci->upper);
|
|
free(ci->lower);
|
|
}
|
|
|
|
double
|
|
pf(q, df1, df2)
|
|
double q, df1, df2;
|
|
{
|
|
double ibeta();
|
|
|
|
return(ibeta(q*df1/(df2+q*df1), df1/2, df2/2));
|
|
}
|
|
|
|
double
|
|
qt(p, df)
|
|
double p, df;
|
|
{
|
|
double t, invibeta();
|
|
|
|
t = invibeta(fabs(2*p-1), 0.5, df/2);
|
|
return((p>0.5?1:-1) * sqrt(t*df/(1-t)));
|
|
}
|
|
|
|
/**********************************************************************/
|
|
/*
|
|
* Incomplete beta function.
|
|
* Reference: Abramowitz and Stegun, 26.5.8.
|
|
* Assumptions: 0 <= x <= 1; a,b > 0.
|
|
*/
|
|
#define DOUBLE_EPS 2.2204460492503131E-16
|
|
#define IBETA_LARGE 1.0e30
|
|
#define IBETA_SMALL 1.0e-30
|
|
|
|
double
|
|
ibeta(x, a, b)
|
|
double x, a, b;
|
|
{
|
|
int flipped = 0, i, k, count;
|
|
double I, temp, pn[6], ak, bk, next, prev, factor, val;
|
|
|
|
if (x <= 0)
|
|
return(0);
|
|
if (x >= 1)
|
|
return(1);
|
|
|
|
/* use ibeta(x,a,b) = 1-ibeta(1-x,b,a) */
|
|
if ((a+b+1)*x > (a+1)) {
|
|
flipped = 1;
|
|
temp = a;
|
|
a = b;
|
|
b = temp;
|
|
x = 1 - x;
|
|
}
|
|
|
|
pn[0] = 0.0;
|
|
pn[2] = pn[3] = pn[1] = 1.0;
|
|
count = 1;
|
|
val = x/(1.0-x);
|
|
bk = 1.0;
|
|
next = 1.0;
|
|
do {
|
|
count++;
|
|
k = count/2;
|
|
prev = next;
|
|
if (count%2 == 0)
|
|
ak = -((a+k-1.0)*(b-k)*val)/
|
|
((a+2.0*k-2.0)*(a+2.0*k-1.0));
|
|
else
|
|
ak = ((a+b+k-1.0)*k*val)/
|
|
((a+2.0*k)*(a+2.0*k-1.0));
|
|
pn[4] = bk*pn[2] + ak*pn[0];
|
|
pn[5] = bk*pn[3] + ak*pn[1];
|
|
next = pn[4] / pn[5];
|
|
for (i=0; i<=3; i++)
|
|
pn[i] = pn[i+2];
|
|
if (fabs(pn[4]) >= IBETA_LARGE)
|
|
for (i=0; i<=3; i++)
|
|
pn[i] /= IBETA_LARGE;
|
|
if (fabs(pn[4]) <= IBETA_SMALL)
|
|
for (i=0; i<=3; i++)
|
|
pn[i] /= IBETA_SMALL;
|
|
} while (fabs(next-prev) > DOUBLE_EPS*prev);
|
|
factor = a*log(x) + (b-1)*log(1-x);
|
|
factor -= gamma(a+1) + gamma(b) - gamma(a+b);
|
|
I = exp(factor) * next;
|
|
return(flipped ? 1-I : I);
|
|
}
|
|
|
|
/*
|
|
* Rational approximation to inverse Gaussian distribution.
|
|
* Absolute error is bounded by 4.5e-4.
|
|
* Reference: Abramowitz and Stegun, page 933.
|
|
* Assumption: 0 < p < 1.
|
|
*/
|
|
|
|
static double num[] = {
|
|
2.515517,
|
|
0.802853,
|
|
0.010328
|
|
};
|
|
|
|
static double den[] = {
|
|
1.000000,
|
|
1.432788,
|
|
0.189269,
|
|
0.001308
|
|
};
|
|
|
|
double
|
|
invigauss_quick(p)
|
|
double p;
|
|
{
|
|
int lower;
|
|
double t, n, d, q;
|
|
|
|
if(p == 0.5)
|
|
return(0);
|
|
lower = p < 0.5;
|
|
p = lower ? p : 1 - p;
|
|
t = sqrt(-2 * log(p));
|
|
n = (num[2]*t + num[1])*t + num[0];
|
|
d = ((den[3]*t + den[2])*t + den[1])*t + den[0];
|
|
q = lower ? n/d - t : t - n/d;
|
|
return(q);
|
|
}
|
|
|
|
/*
|
|
* Inverse incomplete beta function.
|
|
* Assumption: 0 <= p <= 1, a,b > 0.
|
|
*/
|
|
|
|
double
|
|
invibeta(p, a, b)
|
|
double p, a, b;
|
|
{
|
|
int i;
|
|
double ql, qr, qm, qdiff;
|
|
double pl, pr, pm, pdiff;
|
|
double invibeta_quick(), ibeta();
|
|
|
|
/* MEANINGFUL(qm);*/
|
|
qm = 0;
|
|
if(p == 0 || p == 1)
|
|
return(p);
|
|
|
|
/* initialize [ql,qr] containing the root */
|
|
ql = qr = invibeta_quick(p, a, b);
|
|
pl = pr = ibeta(ql, a, b);
|
|
if(pl == p)
|
|
return(ql);
|
|
if(pl < p)
|
|
while(1) {
|
|
qr += 0.05;
|
|
if(qr >= 1) {
|
|
pr = qr = 1;
|
|
break;
|
|
}
|
|
pr = ibeta(qr, a, b);
|
|
if(pr == p)
|
|
return(pr);
|
|
if(pr > p)
|
|
break;
|
|
}
|
|
else
|
|
while(1) {
|
|
ql -= 0.05;
|
|
if(ql <= 0) {
|
|
pl = ql = 0;
|
|
break;
|
|
}
|
|
pl = ibeta(ql, a, b);
|
|
if(pl == p)
|
|
return(pl);
|
|
if(pl < p)
|
|
break;
|
|
}
|
|
|
|
/* a few steps of bisection */
|
|
for(i = 0; i < 5; i++) {
|
|
qm = (ql + qr) / 2;
|
|
pm = ibeta(qm, a, b);
|
|
qdiff = qr - ql;
|
|
pdiff = pm - p;
|
|
if(fabs(qdiff) < DOUBLE_EPS*qm || fabs(pdiff) < DOUBLE_EPS)
|
|
return(qm);
|
|
if(pdiff < 0) {
|
|
ql = qm;
|
|
pl = pm;
|
|
} else {
|
|
qr = qm;
|
|
pr = pm;
|
|
}
|
|
}
|
|
|
|
/* a few steps of secant */
|
|
for(i = 0; i < 40; i++) {
|
|
qm = ql + (p-pl)*(qr-ql)/(pr-pl);
|
|
pm = ibeta(qm, a, b);
|
|
qdiff = qr - ql;
|
|
pdiff = pm - p;
|
|
if(fabs(qdiff) < 2*DOUBLE_EPS*qm || fabs(pdiff) < 2*DOUBLE_EPS)
|
|
return(qm);
|
|
if(pdiff < 0) {
|
|
ql = qm;
|
|
pl = pm;
|
|
} else {
|
|
qr = qm;
|
|
pr = pm;
|
|
}
|
|
}
|
|
|
|
/* no convergence */
|
|
return(qm);
|
|
}
|
|
|
|
/*
|
|
* Quick approximation to inverse incomplete beta function,
|
|
* by matching first two moments with the Gaussian distribution.
|
|
* Assumption: 0 < p < 1, a,b > 0.
|
|
*/
|
|
|
|
static double
|
|
misc_fmin(a, b)
|
|
double a, b;
|
|
{
|
|
return(a < b ? a : b);
|
|
}
|
|
|
|
static double
|
|
misc_fmax(a, b)
|
|
double a, b;
|
|
{
|
|
return(a > b ? a : b);
|
|
}
|
|
|
|
double
|
|
invibeta_quick(p, a, b)
|
|
double p, a, b;
|
|
{
|
|
double x, m, s, invigauss_quick();
|
|
|
|
x = a + b;
|
|
m = a / x;
|
|
s = sqrt((a*b) / (x*x*(x+1)));
|
|
return(misc_fmax(0.0, misc_fmin(1.0, invigauss_quick(p)*s + m)));
|
|
}
|
|
|
|
typedef double doublereal;
|
|
typedef int integer;
|
|
|
|
void
|
|
Recover(a, b)
|
|
char *a;
|
|
int *b;
|
|
{
|
|
printf(a,b);
|
|
exit(1);
|
|
}
|
|
|
|
void
|
|
Warning(a, b)
|
|
char *a;
|
|
int *b;
|
|
{
|
|
printf(a,b);
|
|
}
|
|
|
|
/* d1mach may be replaced by Fortran code:
|
|
mail netlib@netlib.bell-labs.com
|
|
send d1mach from core.
|
|
*/
|
|
|
|
#include <float.h>
|
|
|
|
doublereal F77_SUB(d1mach) (i)
|
|
integer *i;
|
|
{
|
|
switch(*i){
|
|
case 1: return DBL_MIN;
|
|
case 2: return DBL_MAX;
|
|
case 3: return DBL_EPSILON/FLT_RADIX;
|
|
case 4: return DBL_EPSILON;
|
|
case 5: return log10(FLT_RADIX);
|
|
default: Recover("Invalid argument to d1mach()", 0L);
|
|
}
|
|
}
|
|
|