|
|
Boost : |
From: Paul A. Bristow (boost_at_[hidden])
Date: 2001-09-05 08:26:06
> -----Original Message-----
> From: Eric Ford [mailto:eford_at_[hidden]]
> Sent: Tuesday, September 04, 2001 11:16 PM
> To: boost_at_[hidden]
> Subject: [boost] Re: Math Functions
>
> > I am working on some math functions like erf, beta, gamma,
> > and have some questions on preferred practice (hoping to
> > avoid lots of criticisms at a late stage)!
>
> I also recently implemented such functions and had planned to try to
> clean them up, but haven't gotten around to dealing with lots of the
> details, several of which you mentioned in your email. If you're
> interested in seeing what I've got (not documented well),
> http://www.astro.princeton.edu/~eford/math/ has some files.
>
> They include some tempaltes for dealing with inifite sequences,
> infinite series, and continued fractions, which could be used quite
> generally. A also made a template for manipulating numbers by using
> their log that was useful for gamma and beta functions. I've only
> implemented log gamma, gamma, log beta, beta, incomplete beta.
>
I've perused these briefly and declare I am completely out-templated!
My real needs (and those wanting to do some stats) are simpler:
double pNorm01(double z); // Probability of quantile z.
double qNorm01(double p); // Quantile of probability p.
double pT(double t, double df, double ncp = 0.); // Probability of t
double qT(double p, double df, double ncp = 0.); // Quantile of probability
p
double pChisq(double x, double df, double ncp = 0.); // Probability of chisq
double qChisq(double p, double df, double ncp = 0.); // Quantile of
probability p
double pBeta(double x, double a, double b); // Probability of t
double qBeta(double p, double a, double b); // Quantile of probability p
double pF(double f, double dfn, double dfd, double ncp = 0.); // Probability
of t
double qF(double p, double dfn, double dfd, double ncp = 0.); // Quantile of
probability p
double pBinom(double x, double n, double pr); // Probability of t
double pnBinom(double x, double n, double pr); // Quantile of probability p
double pPoisson(double x, double lambda); // Probability of x
double qPoisson(double p, double lambda); // Quantile of probability p
double pGamma(double x, double shape, double scale); // Probability of t
double qGamma(double p, double shape, double scale); // Quantile of
probability p
But of course you need more functions to achieve these, for example
Brown's C version of Morris and Didnato's FORTRAN version (UGH!!!) provides
these:
(Note the automatic translation F2C means that most use double* to pass
values in AND OUT. Seriously un-nice!) I might get round to making them
more C++-ish).
Paul
Dr Paul A Bristow, hetp Chromatography
Prizet Farmhouse
Kendal, Cumbria
LA8 8AB UK
+44 1539 561830
Mobile +44 7714 33 02 04
mailto:pbristow_at_[hidden]
double algdiv(double* a, double* b);
// Computation of ln(gamma(b)/gamma(a+b)) when b >= 8
// in this algorithm, del(x) is the function defined by
// ln(gamma(x)) = (x - 0.5)*ln(x) - x + 0.5*ln(2*pi) + del(x).
double alngam(double* x);
// Returns the natural logarithm of GAMMA(X).
// X --> value at which scaled log gamma is to be returned.
// If X <= 6.0, then use recursion to get X below 3
// then apply rational approximation number 5236, page 246 of
// Hart John F., et al, Computer Approximations, John Wiley, NY,
// SIAM series in Applied Mathematics,
// ISBN 047135601 (1968), reprinted ISBN 0 88275642 7 (1978).
// ISBN 0X04039254 J Wiley (1968)
// If X > 6.0, then use recursion to get X to at least 12 and
// then use formula 5423 of the same source.
double alnrel(double* a);
// Evaluation of the function natural log ln(1 + A)
double apser(double* a,double* b,double* x,double* eps);
// APSER yields the incomplete beta ratio I(SUB(1-X))(B,A) FOR
// A <= MIN(EPS,EPS*B), B*X <= 1, AND X <= 0.5,
// used when A Is Very Small.
double basym(double* a, double* b, double* lambda, double* eps);
// Asymptotic expansion for Ix(a,b) for large a and b.
// lambda = (a + b)*y - b and eps is the tolerance used.
// Assumes that lambda is nonnegative and that
// a and b >= 15.
double bcorr(double* a0,double* b0);
// Evaluation of del(a0) + del(b0) - del(a0 + b0) where
// ln(gamma(a)) = (a - 0.5)*ln(a) - a + 0.5*ln(2*pi) + del(a).
// assumes that a0 >= 8 and b0 >= 8.
double betaln(double* a0, double* b0);
// Evaluation of the logarithm of the BETA function.
double bfrac(double* a,double* b,double* x,double* y,double* lambda, double*
eps);
// Continued fraction expansion for Ix(A,B) when A, B > 1.
// It is assumed that lambda = (A + B)*Y - B.
void bgrat(double* a, double* b, double* x, double* y, double* w, double*
eps, int* ierr);
// Asymptotic expansion for Ix(a,b) when a is larger than b.
// the result of the expansion is added to w.
// Assumes that a >= 15 and b <= 1. eps is the tolerance used.
double bpser(double* a,double* b,double* x,double* eps);
// Power series expansion for evaluating ix(a,b)
// when b <= 1 or b*x <= 0.7. eps is the tolerance used.
void bratio(double* a,double* b,double* x,double* y,double* w, double*
w1,int* ierr);
// Evaluation Of The Incomplete Beta Function Ix(A,B)
// assumed that a and b are nonnegative, and that x <= 1
// and y = 1 - x.
// bratio assigns w and w1 the values
//
// W = Ix(A,B) - return
// W1 = 1 - IX(A,B) - complement.
double brcmp1(int* mu,double* a,double* b,double* x,double* y);
// Evaluation of exp(mu) * (x**a*y**b/beta(a,b))
double brcomp(double* a,double* b,double* x,double* y);
// Evaluation of x**a*y**b/beta(a,b)
double bup(double* a,double* b,double* x,double* y,int* n,double* eps);
// Evaluation of ix(a,b) - ix(a+n,b) where n is a positive integer.
// eps is the tolerance used.
void cdfbet(int* which,double* p,double* q,double* x,double* y,
double* a,double* b,int* status,double* bound);
// Cumulative Distribution Function BETa Distribution.
void cdfbin(int* which,double* p,double* q,double* s,double* xn,
double* pr,double* ompr,int* status,double* bound);
// Cumulative Distribution Function BINomial distribution.
void cdfchi(int* which,double* p,double* q,double* x,double* df,int*
status,double* bound);
// Cumulative Distribution Function - CHI-Square distribution.
void cdfchn(int* which,double* p,double* q,double* x,double* df,
double* pnonc,int* status,double* bound);
// Cumulative Distribution Function - Non-central Chi-Square.
void cdff(int* which,double* p,double* q,double* f,double* dfn,
double* dfd,int* status,double* bound);
// Cumulative Distribution Function - F distribution.
void cdffnc(int* which,double* p,double* q,double* f,double* dfn,
double* dfd,double* phonc,int* status,double* bound);
// Cumulative Distribution Function - Non-central F distribution.
void cdfgam(int* which,double* p,double* q,double* x,double* shape,
double* scale,int* status,double* bound);
// Cumulative Distribution Function - GAMma Distribution.
void cdfnbn(int* which,double* p,double* q,double* s,double* xn,
double* pr,double* ompr,int* status,double* bound);
// Cumulative Distribution Function - Negative BiNomial distribution
void cdfnor(int* which,double* p,double* q,double* x,double* mean,
double* sd,int* status,double* bound);
// Cumulative Distribution Function - NORmal distribution.
void cdfpoi(int* which,double* p,double* q,double* s,double* xlam, int*
status,double* bound);
// Cumulative Distribution Function - POIsson distribution
void cdft(int* which,double* p,double* q,double* t,double* df, int*
status,double* bound);
// Cumulative Distribution Function - T distribution.
void cdftnc(int* which,double* p,double* q,double* t,double* df,
double* pnonc,int* status,double* bound);
// Cumulative Distribution Function - Non-Central T distribution.
void cumbet(double* x,double* y,double* a,double* b,double* cum,double*
ccum);
// Cumulative incomplete BETa distribution.
void cumchi(double* x,double* df,double* cum,double* ccum);
// CUMulative of the CHi-square distribution.
void cumchn(double* x,double* df,double* pnonc,double* cum, double* ccum);
// CUMulative of the Non-central CHi-square distribution.
void cumf(double* f,double* dfn,double* dfd,double* cum,double* ccum);
// CUMulative F distribution.
void cumfnc(double* f,double* dfn,double* dfd,double* pnonc, double*
cum,double* ccum);
// F Non-Central Distribution
// Computes non-central f distribution with dfn and dfd
// Degrees of freedom and noncentrality parameter pnonc.
void cumnbn(double* s,double* xn,double* pr,double* ompr, double*
cum,double* ccum);
// Returns the probability that it there will be S or fewer failures
// before there are XN successes, with each binomial trial having
// a probability of success PR.
void cumnor(double* arg,double* result,double* ccum);
// Computes the cumulative of the normal distribution,
// i.e.,the integral from -infinity to x of
// (1/sqrt(2*pi)) exp(-u*u/2) du
void cumt(double* t,double* df,double* cum,double* ccum);
// CUMulative T-distribution.
void cumtnc(double* t,double* df,double* pnonc,double* cum, double* ccum);
// CUMulative Non-Central T-distribution
// Computes the integral from -infinity to T of the non-central t-density.
double devlpl(double a[], int* n, double* x);
// EVALuate a PoLynomial at X, returns A(1) + A(2)*X + ... + A(N)*X**(N-1)
double dinvnr(double* p, double* q);
// NoRmal distribution INVerse
// Returns X such that CUMNOR(X) = P, i.e., the integral from -
// infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P
void dinvr(int* status,double* x,double* fx, unsigned long *qleft,unsigned
long *qhi);
// Bounds the zero of the function and invokes zror - Reverse Communication
void dstinv(double* zsmall,double* zbig,double* zabsst,
double* zrelst,double* zstpmu,double* zabsto,
double* zrelto);
// SeT INverse finder - (Reverse Communication).
double dt1(double* p,double* q,double* df);
// Approximation toINVerse of the cumulative T distribution.
void dzror(int* status,double* x,double* fx,double* xlo,
double* xhi,unsigned long *qleft,unsigned long *qhi);
// ZeRo of a function -- Reverse Communication
// Performs the zero finding.
// STZROR must have been called before
// this routine in order to set its parameters.
void dstzr(double* zxlo,double* zxhi,double* zabstl,double* zreltl);
// SeT ZeRo finder - Reverse communication version.
// Sets quantities needed by ZROR. The function of ZROR
// and the quantities set is given here.
double erf1(double* x);
// Evaluation of the real error function.
double erfc1(int* ind,double* x);
// Evaluation of the complementary error function
// IF IND = 0 ERFC1(IND, X) = ERFC(X)
// else ERFC1(IND, X) = EXP(X*X)*ERFC(X).
double esum(int* mu,double* x);
// Evaluation of exp(MU + X)
double exparg(int* l);
// largest positive or negative w for which exp(w) can be computed.
double fpser(double* a,double* b,double* x,double* eps);
// Evaluation of ix (a,b) for b < min(eps,eps*a) and x <= 0.5.
double gam1(double* a);
// Computation of 1/gamma(a+1) - 1 for -0.5 <= a <= 1.5
void gaminv(double* a,double* x,double* x0,double* p,double* q, int* ierr);
// Inverse incomplete gamma ratio function
// given positive a, and nonegative p and q where p + q = 1.
// then x is computed where p(a,x) = p and q(a,x) = q.
double gamln(double* a);
// Evaluation of ln(gamma(a)) for positive a.
double gamln1(double* a);
// Evaluation of ln(gamma(1 + a)) for -0.2 <= a <= 1.25
double Xgamm(double* a);
// Evaluation of the gamma function for real arguments.
void grat1(double* a,double* x,double* r,double* p,double* q, double* eps);
// Evaluation of the incomplete gamma ratio functions p(a,x) and q(a,x)
void gratio(double* a,double* x,double* ans,double* qans,int* ind);
// Evaluation of the incomplete gamma ratio functions p(a,x) and q(a,x)
double gsumln(double* a,double* b);
// EVALUATION OF THE FUNCTION LN(GAMMA(A + B))
// FOR 1 <= A <= 2 AND 1 <= B <= 2
double psi(double* xx);
// Evaluation of the digamma function.
double rcomp(double* a,double* x);
// EVALUATION OF EXP(-X)*X**A/GAMMA(A)
double rexp(double* x);
// Evaluation of the function exp(x) - 1
double rlog(double* x);
// Computation of x - 1 - ln(x)
double rlog1(double* x);
// Evaluation of the function x - ln(1 + x)
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk