#include <stdio.h> #include <stdlib.h> #include <math.h> double gammln (double xx); double gammp(double x, double a, double b); double ncbeta(double x, double a, double b, double lam); double betai (double a, double b, double x); double ncf (double f, double n1, double n2, double lam); double probf (double f, double n1, double n2); double probnormi (double u); double probnorm (double x); double ncfinv (double p, double v1, double v2, double lam); double probfi (double p, double n1, double n2); double nct (double t, double df, double nc); double probt (double t, double n1); double probti (double p, double df); double ncti (double p, double df, double nc); double ncgamma (double x, double a, double b, double lam); double probchi (double x, double df); double ncchi (double x, double df, double lam); double ncgammai (double p, double a, double b, double lam); double ncchii (double p, double a, double lam); double probchii (double p, double df);
All of the probability density functions are in Abramowitz and Stegun, so I won't bother with typing them in. They are pretty, but most of the people reading this don't have LATEX, so it will look bad in this document. Besides, leafing through AS is a good adventure. It always reminds me of the times when people had to think about the potential value of their equations before computing them because computation was extremely expensive. These days people just crunch numbers without really understanding the mathematical foundation of the powerful programs we use. Something elegant is being forgotten.
In the following table, I will indicate that the routine comes from Numerical Recipes by NR and the page number. You will need their book if you want to understand the code. If the authors of this reference object to my use of their code, I will remove it from YAMP. AS stands for an Abramowitz and Stegun formula.
Throughout the code, x is for quantiles, p or
u are for probabilities, a or b are
for α or β, nc or lam are for
non–centrality parameters, and df is for degrees of
freedom. The degrees of freedom are continuous throughout.
The following table contains the range of the parameters in the
functions and a short description.
| Function | Description | Ranges |
| gammln(xx) | log Gamma function(NR 168) | xx > 0 |
| gammp(x, a, b) | Incomplete Gamma(NR 172) | x≥0, a, b > 0 |
| ncbeta(u,a,b,lam) | Non-central Beta | u∈[0, 1], a, b > 0, and lam≥0, |
| ncf(x, n1,n2, lam) | Non-central F | ncbeta( u = ax/(ax + b), a = n1/2, b = n2/2) |
| probf(x,n1,n2) | Central F | ncf with lam=0 |
| probnormi( u ) | Inverse Std. Normal (AS 26.2.23) | 0 < u < 1 |
| probnorm( x ) | Standard Normal (NR 175) | x∈ℜ |
| ncfinv( p, v1, v2, lam) | Inverse Non-Central F | p∈(0, 1) |
| nct( t, df, nc) | Non-central t | t∈ℜ, df > 0, nc≥ 0 |
| probt( t, df ) | Central t | same as nct but with nc = 0 |
| probt( p, df ) | Inverse Central t | p∈(0, 1) |
| probti( p, df, nc) | Inverse Non–central t | same as probt but nc > 0 |
| ncgamma(g, a, b, lam) | Non-central Gamma | g, lam≥ 0 and a, b > 0 |
| probchi(x,n1) | Central Chi-Square | ncgamma( g = x/b, a = n1/2, b = 0.5, lam = 0) |
| ncchi(x,n1,lam) | Non-Central Chi-Square | probchi with lam > 0 |
| ncgammai( p,a,b,lam ) | Inverse Non-central gamma | p∈(0, 1) |
| ncchii( p, a, lam ) | Inverse Non-central chi-square | p∈(0, 1) |
| probchii(p, df) | Inverse central Chi-square | p∈(0, 1) |