This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Re: gsl-chisquare cumulative density function



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 = &alpha;

  //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


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]