|
Документ взят из кэша поисковой машины. Адрес
оригинального документа
: http://www.stsci.edu/stsci/meetings/irw/proceedings/snyderd.dir/sectionstar3_5.html
Дата изменения: Mon Apr 18 18:21:02 1994 Дата индексирования: Sun Dec 23 20:55:22 2007 Кодировка: Поисковые слова: п п п п п п п п п п п п п п п п п п п п п п п п п п |
#include <stdio.h>
#include <math.h>
#define EPSILON 0.000001
/* Number of Picard iterations to use */
#define NUM_PICARD 3
float the_function();
double sp_approx(), newton_iteration(), pade(), picard();
/* As listed below, the_function uses the Pade approximation. To use */
/* the Picard iteration, change the "pade" to "picard" in the code */
/* for the_function */
/* Function: the_function */
/* Evaluates F given v, mu, and sigmasquared, where v = r - m */
/* Returns evaluation of F */
/* Remember the normalizing 2*PI factor has been left out */
float the_function(v, mu, sigmasquared)
float v, mu, sigmasquared;
{
double log_q0, log_q1, saddlepoint;
float v_minus_1;
v_minus_1 = v - 1.0;
saddlepoint = pade(v,mu,sigmasquared);
saddlepoint = newton_iteration(v,mu,sigmasquared,saddlepoint);
log_q0 = sp_approx(v,mu,sigmasquared,saddlepoint);
saddlepoint = newton_iteration(v_minus_1,mu,sigmasquared,saddlepoint);
log_q1 = sp_approx(v_minus_1,mu,sigmasquared,saddlepoint);
return((float) (mu * exp(log_q1 - log_q0)));
} /* end the_function */
/* Function: sp_approx */
/* Returns the saddlepoint approximation to the the log of */
/* p(x,mu,sigmasquared) given the saddlepoint found by the */
/* Newton iteration. Remember 2*PI factor has been left out. */
double sp_approx(x,mu,sigmasquared,saddlepoint)
float x, mu, sigmasquared;
double saddlepoint;
{
double mu_exp_minus_s, phi2;
mu_exp_minus_s = mu * exp (-saddlepoint);
phi2 = sigmasquared + mu_exp_minus_s;
return(-mu + (saddlepoint * (x + 0.5 * sigmasquared * saddlepoint))
+ mu_exp_minus_s - 0.5 * log(phi2));
} /* end sp_approx */
/* Function: newton_iteration */
/* Returns the saddlepoint found by Newton iteration for a given */
/* x, mu, sigmasquared and an initial esimate of the saddlepoint */
/* (found with either the Pade or Picard approaches */
double newton_iteration(x, mu, sigmasquared, initial_saddlepoint)
float x, mu, sigmasquared;
double initial_saddlepoint;
{
double mu_exp_minus_s, saddlepoint, change;
saddlepoint = initial_saddlepoint;
do {
mu_exp_minus_s = mu * exp(-saddlepoint);
change = (x + sigmasquared * saddlepoint - mu_exp_minus_s)
/ (sigmasquared + mu_exp_minus_s);
saddlepoint -= change;
} while((fabs(change) > EPSILON * fabs(saddlepoint)));
return(saddlepoint);
} /* end newton_iteration */
/* Function: pade */
/* Returns the initial saddlepoint estimated by the Pade approx. */
double pade(x, mu, sigmasquared)
float x, mu, sigmasquared;
{
double bterm;
bterm = x - 2*sigmasquared - mu;
return(-log(0.5 * (bterm
+ sqrt(bterm*bterm + 4*mu*(2*sigmasquared + x))) / mu));
} /* end pade */
/* Function: picard */
/* Returns the initial saddlepoint estimated by the Picard iter. */
double picard(x, mu, sigmasquared)
float x, mu, sigmasquared;
{
int i;
double argument_to_log, saddlepoint, taylor;
/* Use Taylor approx. to get starting point for Picard iteration. */
saddlepoint = taylor = (mu - x) / (mu + sigmasquared);
for (i = 0; i < NUM_PICARD; i++)
{
argument_to_log = mu / (x + sigmasquared * saddlepoint);
if (argument_to_log <= 0.0) /* Break out of loop if */
return(taylor); /* arg. to log goes neg.*/
else saddlepoint = log(argument_to_log);
}
return(saddlepoint);
} /* end picard */
Acknowledgments
This work was supported in part by the National Science Foundation under grant number MIP-9101991.