This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: gsl-chisquare cumulative density function
- To: bjg at network-theory dot co dot uk
- Subject: Re: gsl-chisquare cumulative density function
- From: "Varghese John" <vjohn111 at hotmail dot com>
- Date: Thu, 27 Sep 2001 08:28:52 -0400
- Cc: gsl-discuss at sources dot redhat dot com
Hi,
Not sure if this will be of help, but here are three solutions.
1)the Normal CDF is nothing but
the error function upto normalisation.
So you can use the gsl Error function : gsl_sf_erfc.
2) Use a simple rational approx.
//------------------------------------------------------------
//Rational Approx to NormalCDF-- Good to four decimals(?)
//------------------------------------------------------------
float normCdf(float x){
float a1=.319381530;
float a2= -0.356563782;
float a3=1.781477937;
float a4=-1.821255978;
float a5=1.330274429;
float g=0.2316419;
float k=1/(1+g*x);
float m=1/(1-g*x);
if(x >= 0 ){
return 1-(1/sqrt(2*3.141))*exp(-x*x/2)*(k*a1+a2*pow(k,2) + a3*pow(k,3) +
a4*pow(k,4) +a5*pow(k, 5));
}
else{
return (1/sqrt(2*3.141))*exp(-x*x/2)*(m*a1+a2*pow(m,2) + a3*pow(m,3) +
a4*pow(m,4) +a5*pow(m, 5));
}
}
//------------------------------------------------------------
3)Use the numerical Integration.
I also have tried using the Numerical Integration routine
to get good results.
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_integration.h>
#include <math.h>
double f (double x, void * params) {
double alpha = *(double *) params ;
double f= (1/(M_SQRT2*M_SQRTPI))*exp(-x*x/2);
return f ;
}
int main ()
{
//Numerical integration
gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);
double result, error;
double expected = 0.1587;
double alpha = 1.0;
double a = 2;
gsl_function F;
F.function = &f;
F.params = α
//gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000, w, &result, &error);
//integrate pdf from a to infinity.
gsl_integration_qagiu (&F, a, 0, 1e-7, 1000, w, &result, &error);
printf("NORMCDF(2) = % .18f\n", 1-result);
//printf("result = % .18f\n", result);
//printf("exact result = % .18f\n", expected);
//printf("estimated error = % .18f\n", error);
//printf("actual error = % .18f\n", result - expected);
//printf("intervals = %d\n", w->size);
return(0);
}
_________________________________________________________________
Get your FREE download of MSN Explorer at http://explorer.msn.com/intl.asp