Re: [Rd] pgamma discontinuity (PR#7307)

From: <terra_at_gnome.org>
Date: Wed 27 Oct 2004 - 16:15:38 EST

For the record, this is what I am going to use in Gnumeric for now.

Morten

/*
/*
     do {
	  a += 1.;
	  c *= x / a;
	  sum += c;
     } while (c > DBL_EPSILON * sum);

     return sum;

}
/*
     a = 1. - alph;
     b = a + x + 1.;
     pn1 = 1.;
     pn2 = x;
     pn3 = x + 1.;
     pn4 = x * b;
     sum = pn3 / pn4;
     for (n = 1; ; n++) {
	  a += 1.;/* =   n+1 -alph */
	  b += 2.;/* = 2(n+1)-alph+x */
	  an = a * n;
	  pn5 = b * pn3 - an * pn1;
	  pn6 = b * pn4 - an * pn2;
	  if (fabs(pn6) > 0.) {
	       osum = sum;
	       sum = pn5 / pn6;
	       if (fabs(osum - sum) <= DBL_EPSILON * fmin2(1., sum))
		    break;
	  }
	  pn1 = pn3;
	  pn2 = pn4;
	  pn3 = pn5;
	  pn4 = pn6;
	  if (fabs(pn5) >= xlarge) {
	       /* re-scale the terms in continued fraction if they are large */
#ifdef DEBUG_p
	       REprintf(" [r] ");
#endif
	       pn1 /= xlarge;
	       pn2 /= xlarge;
	       pn3 /= xlarge;
	       pn4 /= xlarge;
	  }
     }

     return sum;

}

double pgamma(double x, double alph, double scale, int lower_tail, int log_p) {

     double arg;

#ifdef IEEE_754

     if (ISNAN(x) || ISNAN(alph) || ISNAN(scale))
	  return x + alph + scale;
#endif
     if(alph <= 0. || scale <= 0.)
	  ML_ERR_return_NAN;
     if (x <= 0.)
	  return R_DT_0;
     x /= scale;
#ifdef IEEE_754
     if (ISNAN(x)) /* eg. original x = scale = +Inf */
	  return x;

#endif
     if (x < 0.99 * (alph + 1000) && x < 10 + alph && x < 10 * alph) {
	  /*
	   * The first condition makes sure that terms in the sum get at
	   * least 1% smaller, no later than after 1000 terms.
	   *
	   * The second condition makes sure that the terms start getting
	   * smaller (even if just a tiny bit) after no more than 10 terms.
	   *
	   * The third condition makes sure that the terms don't get huge
	   * before they start getting smaller.
	   */
	  double C1mls2p = /* 1 - log (sqrt (2 * M_PI)) */
	       0.081061466795327258219670263594382360138602526362216587182846;
	  double logsum = log1p (pgamma_series_sum (x, alph));
	  /*
	   * The first two terms here would cause cancellation for extremely
	   * large and near-equal x and alph.  The first condition above
	   * prevents that.
	   */
	  arg = (alph - x) + alph * log (x / (alph + 1)) + C1mls2p -
	       log1p (alph) / 2 - stirlerr (alph + 1) + logsum; 
     } else if (alph < x && alph - 100 < 0.99 * x) {
	  /*
	   * The first condition guarantees that we will start making
	   * a tiny bit of progress right now.
	   *
	   * The second condition guarantees that we will make decent
	   * progress after no more than 100 loops.
	   */
	  double lfrac = log (pgamma_cont_frac (x, alph));
	  arg = lfrac + alph * log(x) - x - lgammafn(alph);
	  lower_tail = !lower_tail;
     } else {
	  /*
	   * Approximation using pnorm.  We only get here for fairly large
	   * x and alph that are within 1% of each other.
	   */
	  double r3m1 = expm1 (log (x / alph) / 3);
	  return pnorm (sqrt (alph) * 3 * (r3m1 + 1 / (9 * alph)),
			0, 1, lower_tail, log_p);
     }

     if (lower_tail)
	  return log_p ? arg : exp(arg);
     else
	  return log_p ? R_Log1_Exp(arg) : -expm1(arg);
}

R-devel@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Wed Oct 27 16:23:54 2004

This archive was generated by hypermail 2.1.8 : Wed 03 Nov 2004 - 22:45:25 EST