R-beta: problem with locfit

[guido@sirio.stat.unipd.it]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._