Hi<br /><br />I am trying to implement a statistical function, similar to the <b>chidist</b> function in MS Excel or OpenOffice. This function uses both the gamma function and the incomplete gamma function. However, when I try to calculate it, my program goes into an endless loop. Here's the code:<br /><br /><!--c1--><div class='codetop'>CODE</div><div class='codemain'><!--ec1--><br />#include <stdio.h><br />#include <math.h><br />#include "StatisticalFunctions.h"<br />#include "HypothesisTesting.h"<br /><br />/* Complement of the error function approximation<br /><br />    erfc(z) = t * exp(-z*z)    * exp( a0 + a1*t + a2*t^2 + a3*t^3 ...<br />  where t = 1/(1 + z/2)<br /><br />  Error:  1.2e-7 maximal relative error for z >= 0<br />  Source:  Numerical Recipies, Chapter 6, Special functions  */<br /><br />double erfc (double z) {<br /><br />    const double a0 = -1.26551223;<br />    const double a1 =  1.00002368;<br />    const double a2 =  0.37409196;<br />    const double a3 =  0.09678418;<br />    const double a4 = -0.18628806;<br />    const double a5 =  0.27886807;<br />    const double a6 = -1.13520398;<br />    const double a7 =  1.48851587;<br />    const double a8 = -0.82215223;<br />    const double a9 =  0.17087277;<br />    const double t = 2./(2. + z);<br />    return t*exp(-z*z+a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*a9))))))))); <br /><br />};<br /><br /><br /><br /><br />/* Gamma function approximation<br /><br />    gamma(z) = ( (sqrt(2*pi)/z) * ( p0 + (p1/z+1) + ... + (p6/z+6) )) * <br />       ( (z+5.5)^(z+0.5) ) * exp( -1*(z+5.5) )<br />    <br />  Error:  Unknown<br />  Source:  Unknown  */<br /><br />double gamma (double z) {<br /><br />    const double p0 = 1.000000000190015;<br />    const double p1 = 76.18009172947146;<br />    const double p2 = -86.50532032941677;<br />    const double p3 = 24.01409824083091;<br />    const double p4 =  -1.231739572450155;<br />    const double p5 = 1.208650973866179*0.001;<br />    const double p6 = -5.395239384953*0.000001;<br /><br />    double P = p0 + (p1/(z+1)) + (p2/(z+2)) + (p3/(z+3)) + (p4/(z+4)) + (p5/(z+5)) + (p6/(z+6));<br />    return (sqrt(2*M_PI)/z)*P*pow(z+5.5,z+0.5)*exp( -1*(z+5.5) );<br /><br />};<br /><br /><br /><br /><br />/* Incomplete gamma function approximation<br /><br />    incgamma(x,a) = (x^a)*exp(-x)*sum(n=0:infinity){(x^n)/(a*(a+1)*...*(a+n))}<br /><br />  Error:  Arbitrary.  As accurate as double allows<br />  Source:  Unknown  */<br /><br />double incgamma (double x, double a) {<br /><br />    printf("Calculating incgamma\n");<br />    double sum = 0;<br />    double term = 1.0/a;<br />    int n = 1;<br />    while (term != 0) {<br />  sum = sum + term;<br />  term = term*(x/(a+n));<br />  n++;<br />    }<br />    return pow(x,a)*exp(-1*x)*sum;<br /><br />}<br /><br /><br /><br />/* Chi-square distribution approximation<br /><br />    chidist(x,v) = ( exp((-1*x)/2) * x^( (v/2) - 1) ) / ( 2^(v/2) * gamma(v/2) )<br />  where x >= 0<br /><br />  Error:  Dependant on the gamma and incomplete gamma functions<br />  Source:  Standard Statistical Function  */<br /><br />double chidist (double x, double v) {<br />    return incgamma(x/2.0,v/2.0)/gamma(x/2.0);<br /><br />}<br /><!--c2--></div><!--ec2--><br /><br />As you can see, I also implemented the error function (erfc). I'm not sure if that is correct either, but I found the source code (in C) somewhere on the web. The algorithms for the gamma and incomplete gamma function I found at <a href="http://www.rskey.org/gamma.htm" target="_blank">Calculators and the Gamma Function</a>, but I had to implement that myself, so I'm sure that is where the problem is. I attempted to rewrite the incomplete gamma function in an iterative manner, but when I call the function the value of the variable 'term' goes to infinity and stays there (so it's obvioustly not going to converge). I'm a little stumped here so any help would be greatly appreciated. The incgamma function is supposed to converge in close to 100 iterations.<br /><br />Thanks in advance
</div>
</div>