Документ взят из кэша поисковой машины. Адрес
оригинального документа
: 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.