gsl-chisquare cumulative density function

[Varghese John]

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