From mboxrd@z Thu Jan  1 00:00:00 1970
From: "Varghese John" <vjohn111@hotmail.com>
To: bjg@network-theory.co.uk
Cc: gsl-discuss@sources.redhat.com
Subject: Re: gsl-chisquare cumulative density function 
Date: Wed, 19 Dec 2001 13:20:00 -0000
Message-id: <F188NuhsZnOXOweys1S0000b058@hotmail.com>
X-SW-Source: 2001/msg00528.html

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