/* Prototypes for various useful statistical subroutines */ /* Version date: April 2, 2007 */ #ifndef STATLIB_HX /* Make sure this is not read more than once */ /* If STATLIB_HX is defined, this causes a jump */ /* to the #endif at the end of this file */ #define STATLIB_HX /* Reading a header file twice can sometimes */ /* cause errors with fussy compilers */ #include #include #include #include #include /* Basic math functions like sqrt() and log() */ #include /* For set_newsignals() */ #include /* For initializing a random-number generator */ extern const char *statlib_version_date; /* Statlib version date in a string */ /* Handy general constants */ #define TRUE 1 #define FALSE 0 #define RET_OK 0 #define RET_ERROR (-1) /* Subroutines to calculate P-values */ /* Gaussian integral: */ /* Returns Prob(N(0,1) <= x) */ double normint(double x); /* Two-sided normal P-value */ /* Returns Prob(|N(0,1)| >= |x|) = 2Prob(N(0,1) >= |x|) = 2*(1-normint(fabs(x))); */ double normpval(double x); /* F-statistics */ /* Returns Prob(F(df1,df2)<=xx) */ /* df1 and df2 are the two degrees of freedom for an F distribution */ double fdist (double df1, double df2, double xx); /* P-value for F distribution */ /* Returns Prob(F(df1,df2)>=xx) */ /* df1 and df2 are the two degrees of freedom for an F distribution */ double pvfstat(double df1, double df2, double xx); /* Quantiles of the F distribution: */ /* Returns xx in Prob(F(df1,df2)<=xx) = pp */ /* Must have df1>0, df2>0, and 0 < pp < 1 */ double invfdist(double df1, double df2, double pp); /* F-distribution critical value: */ /* Returns xx in Prob(F(df1,df2)>=xx) = pp */ /* Must have 0 < pp < 1 */ /* df1 and df2 are the two degrees of freedom for an F distribution */ double crit_fdist(double df1, double df2, double pp); /* T-statistics */ /* Returns Prob( t(df) <= xx) for -infty < xx < +infty */ /* df is the degrees of freedom for a t distribution */ double tdist(double df, double xx); /* Two-sided P-value for t distribution */ /* Returns Prob( |t(df)|>=|x|) = 2Prob(t(df) >= |x|) */ /* df is the degrees of freedom for a t distribution */ double pvttest(double df, double xx); /* Returns xx>=0 in Prob(|t(df)|>=xx) = 2Prob(t(df)>=xx) = pp */ /* Must have 0 < pp < 1 */ /* df is the degrees of freedom for a t distribution */ double crit_tdist(double df, double pp); /* Binomial probabilities */ /* Returns Prob(B(n,p) >= k) */ double binomtail(int k, int n, double p); /* Chi-square and Poisson statistics: */ /* Chi-square P-value */ /* Returns Prob( Chi(df)>=x ) */ /* df is the degrees of freedom for a chi-square distribution */ double pvchisq(double df, double xx); /* Returns xx in Prob(Chi(df)>=xx) = pp */ /* Must have 0 < pp < 1 */ /* df is the degrees of freedom for a chi-square distribution */ double crit_chisq(double df, double pp); /* Returns xx in Prob(Chi(df)<=xx) = pp */ /* Must have 0 < pp < 1 */ /* df is the degrees of freedom for a chi-square distribution */ double invchisq(double df, double pp); /* Returns xx in Prob(Gamma(alf,beta)<=xx) = pp */ /* Must have alf>0, beta>0, and 0 < pp < 1 */ double invgamma(double alpha, double beta, double pp); /* Poisson tail probabilities: */ /* Return Pr(Poi(A) >= n) = Sum(j=n,infty) Pr(Poi(A)=j) */ double poissontail(int n, double x); /* Some utilities that are occasionally useful: */ /* log(Gamma(z)) where Gamma(z) is the gamma function */ double lngam(double z); /* log(Beta(a,b)) where B(a,b) is the beta function */ double lnbeta(double aa, double bb); /* Avoid Microsoft apology windows for certain program errors */ /* (Technically: Provide better command-line behavior for */ /* errors that generate SIGSEGV, SIGFPE, and SIGILL signals) */ void set_newsignals(void); /* Random-number generator (RNG) package */ /* Using a custom RNG means identical (and reliable) output with */ /* the same starting seeds, even if the program is run at different times */ /* on different UNIX, PC, or Mac machines. */ /* Initialize the RNG with seed=uu, or from the system clock if uu=0 */ /* Returns the value of the seed */ /* (This will be different from uu if uu=0). */ unsigned long int setrand(unsigned long int uu); /* Return the current value of the seed */ unsigned long int dispseed(void); /* Display how many random numbers have been used so far */ void show_randnums_used(void); /* Applications of random-number package: */ /* Return a random integer in 0,1,2,...,n-1 */ int nrand(int n); /* Return a random double in (0,1) */ double drand(void); /* `Perfect shuffle' of an array int deck[] of length num */ void shuffle(int *deck, int num); /* `Perfect shuffle' of an array double deck[] of length num */ void shuffle_dub(double *deck, int num); /* `Perfect shuffle' of an array char deck[] of length num */ void shuffle_char(char *deck, int num); /* Return a random variable with an exponential distribution: */ double exprand(void); /* Return a random variable with a standard normal N(0,1) distribution: */ double normrand(void); /* Other useful utilities: */ /* Read an array xx[] of length nn and stores the sample mean */ /* and sample variance in doubles whose ADDRESSES */ /* are in `meanptr' and `varptr' */ /* Set xmeanptr=NULL or xvarptr=NULL if you don't want that */ /* particular value. */ void getmeanvar (const double *xx, int mm, double *xmeanptr, double *xvarptr); double getmedian_sorted (const double *xx, int mm); /* Find midranks for an array of doubles */ /* Assume an array wwdat[] of length ntot */ /* Put midrank of wwdat[i] at mrank[i] */ /* Put tiegroup sizes (>1) in tgsize[] with number at *ntgptr */ void makeranks(double *mrank, const double *wwdat, int ntot, int *tgsize, int *ntgptr); /* Return rank-sum tiesum from (tgsize,ntiegroups) array */ double cubic_tiesum(const int *tgsize, int ntiegroups); /* Provide footnotes for a table */ char *starp(double pp); /* Comparison function for sorting doubles using system `quick sort' */ /* If zarr[] is an array of nn doubles, then */ /* */ /* qsort(zarr,nn, sizeof(double),dubsort); */ /* */ /* will sort zarr[] in increasing order. If iarr[] is an array of nn */ /* integers, */ /* */ /* qsort(iarr,nn, sizeof(int),intsort); */ /* */ /* will sort it to increasing order. */ int sortdub(const void *p, const void *q); int dubsort(const void *p, const void *q); int sortint(const void *p, const void *q); int intsort(const void *p, const void *q); /* Sometimes useful for squares and cubes of complicated expressions */ double square(double x); double cube(double x); /* Show integers in a prettier form */ /* 1,024,076,332 instead of 1024076332 */ /* commnum: 200 2000 20,000 200,000 */ /* commnumf: 200 2,000 20,000 200,000 */ /* Warning: Output is in a character buffer, and is NOT an integer */ char *commnum(int nn); char *commnum_unsigned(unsigned long nn); char *commnumf(int nn); /* Alternatives if several are to be used together */ char *commnum1(int nn); char *commnum2(int nn); #endif /* STATLIB_HX */