Re: gsl-chisquare cumulative density function

Re: gsl-chisquare cumulative density function

[Brian Gough]

Varghese John writes:
 > 
 > 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.
 > 
 > 3)Use the numerical Integration.

Right, ideally I would like to see a separate cdf/ directory with
functions for direct computation of both P(x')=integ(p(x),-inf,x') and
Q(x)=1-P(x) to full precision.  A lot of the cases are already
available in specfunc.

Re: 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 = &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

Re: gsl-chisquare cumulative density function

[Brian Gough]

Gilles Frising writes:
 >  I wanted to use the cumulative chisquare pdf of gsl:
 > `gsl_ran_chisq_pdf'.  I checked this function in maple,cernlib and
 > gsl, and i found different result for gsl.
 >
 >  MAPLE:
 >  statevalf[cdf,chisquare[6]](3.);
 >       .1911531695
 >
 >  CERNLIB:
 >  prob: 0.191153
 >
 >  GSL:
 >  gsl_ran_chisq_pdf(3.,6.): 0.125511
 >                            =========
 >
 >  Is this a bug or have i done sth. wrong?

Hi,

The pdf is the probability distribution function (point value), so
this is a different quantity from the cdf. 

The cumulative cdf's are not implemented in the present version of
GSL.  We are looking for a volunteer to write them.

In the meantime you'd have to look up the formula and use the
corresponding routines from the specfunc directory -- it's an
incomplete gamma function, I believe.

regards
Brian Gough

gsl-chisquare cumulative density function

[Gilles Frising]

Hello there,

I wanted to use the cumulative chisquare pdf of gsl:
`gsl_ran_chisq_pdf'.
I checked this function in maple,cernlib and gsl, and i found different
result for gsl.



MAPLE:

statevalf[cdf,chisquare[6]](3.);
      .1911531695 
       ==========  



CERNLIB:

function: http://wwwinfo.cern.ch/asdoc/shortwrupsdir/g100/top.html

  extern "C" float prob_(float&,int&);
  int n=6;
  float x=3.;
  cout<<"prob: "<<1.-prob_(x,n)<<endl;

prob: 0.191153 
      ========  




GSL:

gsl_ran_chisq_pdf(3.,6.):  0.125511
                           =========


Is this a bug or have i done sth. wrong?

ciao,

Gilles