I installed the locfit package under Linux (gcc 2.7.2). Installation
was ok but> x <- runif(200)
> y.compl <- 10*x*x*rgamma(200,3)
> med.y <- median(y.compl)
> cens <- ifelse(y.compl<=med.y,1,0)
> y <- cens * y.compl + (1-cens)*med.y
> library(locfit)
> m <- locfit(y~x,cens=cens,family="gamma")
/usr/local/src/R-0.61.1/bin/R.binary: can't resolve symbol 'igamma'
I took a look to the original code, not the R/S one.
There igamma (and also ibeta which the r distribution misses too)
are defined in the file random.c.
I inserted them in the locfit/src-c/random.c (see below) and then
things seem to go well (= the previous example works and
give sensible results).
Perhaps, a better possibility is to define igamma and ibeta using
the R pgamma and pbeta function. However, I have not tried it.
guido masarotto
--------------------------------------------------------------------
#define IBETA_LARGE 1.0e30
#define IBETA_SMALL 1.0e-30
#define IGAMMA_LARGE 1.0e30
#define DOUBLE_EPS 2.2204460492503131E-16
double ibeta(x, a, b)
double x, a, b;
{ int flipped = 0, i, k, count;
double I = 0, 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,z) */
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 -= LGAMMA(a+1) + LGAMMA(b) - LGAMMA(a+b);
I = exp(factor) * next;
return(flipped ? 1-I : I);
}
/*
* Incomplete gamma function.
* Reference: Abramowitz and Stegun.
* Assumptions: x >= 0; df > 0.
*/
double igamma(x, df)
double x, df;
{ double factor, term, gintegral, pn[6], rn, ak, bk;
double increment, df1;
int i, count, k;
if (x <= 0.0) return(0.0);
if (df < 1.0)
{ increment = exp(df*log(x) - x - LGAMMA(df + 1.0));
df1 = df + 1.0;
} else
{ increment = 0.0;
df1 = df;
}
factor = exp(df1*log(x) - x - LGAMMA(df1));
if (x > 1.0 && x >= df1)
{ pn[0] = 0.0;
pn[2] = pn[1] = 1.0;
pn[3] = x;
count = 1;
rn = 1.0 / x;
do
{ count++;
k = count / 2;
gintegral = rn;
if (count%2 == 0)
{ bk = 1.0;
ak = (double)k - df1;
} else
{ bk = x;
ak = (double)k;
}
pn[4] = bk*pn[2] + ak*pn[0];
pn[5] = bk*pn[3] + ak*pn[1];
rn = pn[4] / pn[5];
for (i=0; i<4; i++)
pn[i] = pn[i+2];
if (pn[4] > IGAMMA_LARGE)
for (i=0; i<4; i++)
pn[i] /= IGAMMA_LARGE;
} while (fabs(gintegral-rn) > DOUBLE_EPS*rn);
gintegral = 1.0 - factor*rn;
} else
{ gintegral = term = 1.0;
rn = df1;
do
{ rn += 1.0;
term *= x/rn;
gintegral += term;
} while (term > DOUBLE_EPS*gintegral);
gintegral *= factor/df1;
}
return(increment + gintegral);
}
-----------------------------------------------------------------------------
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._