public inbox for gsl-discuss@sourceware.org
 help / color / mirror / Atom feed
* gsl-chisquare cumulative density function
@ 2001-12-19 13:20 Gilles Frising
  2001-12-19 13:20 ` Brian Gough
  0 siblings, 1 reply; 4+ messages in thread
From: Gilles Frising @ 2001-12-19 13:20 UTC (permalink / raw)
  To: gsl-discuss

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

^ permalink raw reply	[flat|nested] 4+ messages in thread
* Re: gsl-chisquare cumulative density function
@ 2001-12-19 13:20 Varghese John
  2001-12-19 13:20 ` Brian Gough
  0 siblings, 1 reply; 4+ messages in thread
From: Varghese John @ 2001-12-19 13:20 UTC (permalink / raw)
  To: bjg; +Cc: gsl-discuss

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

^ permalink raw reply	[flat|nested] 4+ messages in thread

end of thread, other threads:[~2001-12-19 13:20 UTC | newest]

Thread overview: 4+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2001-12-19 13:20 gsl-chisquare cumulative density function Gilles Frising
2001-12-19 13:20 ` Brian Gough
2001-12-19 13:20 Varghese John
2001-12-19 13:20 ` Brian Gough

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).