Chi squared cdf

[Labdabeta]

I need to calculate chi squared cdf from x->infinity with n=255 degrees of freedom. I cannot seem to get it to work. You can find formulae for chi squared distributions using google. Here is my non-functional approximation code:

double riemannsum(double(*fnc)(double),double dx, double xmin, double xmax)
{
    double ret=0;
    double x=xmin;
    while (x<xmax)
    {
        ret+=fnc(x);
        x+=dx;
    }
    return ret;
}
double lowergammaintegrand(double x)
{
    return pow(x,254)*pow(M_E,-x);//I am only using 255 degrees of freedom in this program
}
double lowergamma(double x)
{
    return riemannsum(&lowergammaintegrand,10,x,10000000);
}
double chicdf(double x)
{
    static double gamma1275=3.40511e214;//this is Gamma(127.5)=Gamma(k/2)
    return lowergamma(x)/gamma1275;
}

[histrungalot]

Just curious, why are you integrating over the CDF and not the PDF?
And wouldn't the value of x->infinity = 1 - chicdf(x) (chicdf not being the integral but the value of chicdf at x)?
Just thinking, I could be totally wrong.

[Labdabeta]

I tried integrating over the pdf, but the values were all so small that they were being stored as 0.0, so I used wikipedia to find the integral above. I have been playing around on my TI calculator and have found an approximate expression for chicdf(X):

double chicdf(double x)
{
    return (-1/(1+pow(M_E,(-(x-256)/14))))+1;
}

Its maximum error at x=255 is about -0.3 which is about a 6% error. Any ideas on a better solution?

[histrungalot]

I have never used it but maybe GNU Scientific Library,

http://www.gnu.org/software/gsl/manual/html_node/The-Chi_002dsquared-Distribution.html
or
http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/dist/dist_ref/dists/chi_squared_dist.html

[StuXYZ]

Your errors are dominated by overflow. These kind of functions are normally computed in a logorithmic maniford, not the real number manifold.

Consider: pow(x,254) almost always in overflow for most of your x range.
pow(M_E,-x) at underflow for most of your x range

Those two can be combined: exp( log(x)*256-x), that reduces the over/under flow, and obviously doing the whole problem in the log space reduces it still further.

Your range from x to 100000 is going to do that unless you convert the integration to a log basis, if you calculated log_chi_cdf, you should be ok -- then just convert at the end

[Labdabeta]

Thank you.