Statistics/Chi-squared distribution

Statistics/Chi-squared distribution is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The probability density function (pdf) of the chi-squared (or χ2) distribution as used in statistics is

$f(x;\,k)={\dfrac {x^{{\frac {k}{2}}-1}e^{-{\frac {x}{2}}}}{2^{\frac {k}{2}}\Gamma \left({\frac {k}{2}}\right)}}$, where ${\textstyle x>0}$

Here, ${\textstyle \Gamma (k/2)}$ denotes the Gamma_function.

The use of the gamma function in the equation below reflects the chi-squared distribution's origin as a special case of the gamma distribution.

In probability theory and statistics, the gamma distribution is a two-parameter family of continuous probability distributions. The exponential distribution, Erlang distribution, and chi-square distribution are special cases of the gamma distribution.

The probability density function (pdf) of the gamma distribution is given by the formula

$f(x;k,\theta )={\frac {x^{k-1}e^{-x/\theta }}{\theta ^{k}\Gamma (k)}}\quad {\text{ for }}x>0{\text{ and }}k,\theta >0.$

where Γ(k) is the Gamma_function, with shape parameter k and a scale parameter θ.

The cumulative probability distribution of the gamma distribution is the area under the curve of the distribution, which indicates the increasing probability of the x value of a single random point within the gamma distribution being less than or equal to the x value of the cumulative probability distribution. The gamma cumulative probability distribution function can be calculated as

$F(x;k,\theta )=\int _{0}^{x}f(u;k,\theta )\,du={\frac {\gamma \left(k,{\frac {x}{\theta }}\right)}{\Gamma (k)}},$

where $\gamma \left(k,{\frac {x}{\theta }}\right)$ is the lower incomplete gamma function.

The lower incomplete gamma function can be calculated as

$$x^{s}\,\Gamma (s)\,e^{-x}\sum _{m=0}^{\infty }{\frac {x^{m}}{\Gamma (s+m+1)}}$$

and so, for the chi-squared cumulative probability distribution with chi-squared k, we have, substituting chi-square k into s as k/2,

$$F(x;\,k)={\frac {1}{\Gamma (k)}}x^{(k/2)}\,\Gamma (k/2)\,e^{-x}\sum _{m=0}^{\infty }{\frac {x^{m}}{\Gamma ({\frac {k}{2}}+m+1)}}.$$

In practice, this series formula is often subject to accumulated errors from rounding in the frequently used region where x and k are under 10 and near one another. You may therefore instead use a mathematics function library, if available for your programming task, to calculate gamma and incomplete gamma.

Task

• Calculate and show the values of the χ2(x; k) for k = 1 through 5 inclusive and x integer from 0 and through 10 inclusive.

• Create a function to calculate the cumulative probability function for the χ2 distribution. This will need to be reasonably accurate (at least 6 digit accuracy) for k = 3.

• Calculate and show the p values of statistical samples which result in a χ2(k = 3) value of 1, 2, 4, 8, 16, and 32. (Statistical p values can be calculated for the purpose of this task as approximately 1 - P(x), with P(x) the cumulative probability function at x for χ2.)

The following is a chart for 4 airports:

Flight Delays

Airport	On Time	Delayed	Totals
Dallas/Fort Worth	77	23	100
Honolulu	88	12	100
LaGuardia	79	21	100
Orlando	81	19	100
All Totals	325	75	400
Expected per 100	81.25	18.75	100

χ2 on a 2D table is calculated as the sum of the squares of the differences from expected divided by the expected numbers for each entry on the table. The k for the chi-squared distribution is to be calculated as df, the degrees of freedom for the table, which is a 2D parameter, (count of airports - 1) * (count of measures per airport - 1), which here is (4 - 1 )(2 - 1) = 3.

• Calculate the Chi-squared statistic for the above table and find its p value using your function for the cumulative probability for χ2 with k = 3 from the previous task.

Stretch task

• Show how you could make a plot of the curves for the probability distribution function χ2(x; k) for k = 0, 1, 2, and 3.

Phix

This example is incomplete. just a translation of code on the talk page Please ensure that it meets all task requirements and remove this message.

with javascript_semantics
constant p = {   0.99999999999980993, 
               676.5203681218851,   
             -1259.1392167224028, 
               771.32342877765313,  
              -176.61502916214059,  
                12.507343278686905, 
                -0.13857109526572012,
                 9.9843695780195716e-6,
                 1.5056327351493116e-7 }

function gamma(atom z)
    if z<0.5 then
        return PI / (sin(PI*z)*gamma(1-z))
    end if
    z -= 1;
    atom x := p[1];
    for i=1 to length(p)-1 do x += p[i+1]/(z+i) end for
    atom t = z + length(p) - 1.5;
    return sqrt(2*PI) * power(t,z+0.5) * exp(-t) * x
end function

function pdf(atom x, k)
    if (x <= 0) then return 0 end if // this should be in task description
    return exp(-x/2) * power(x,k/2-1) / (power(2,k/2) * gamma(k/2))
end function

function cpdf(atom x, k)
    x = x / 2  // need to do this to agree with Wikipedia formula for k = 2
    atom t = exp(-x) * power(x,k/2) * gamma(k/2) / gamma(k),
         s = 0,
         m = 0,
         tol = 1e-15 // say
    while true do
        atom term = power(x,m) / gamma(k/2 + m + 1)
        s = s + term
        if abs(term) < tol then exit end if
        m = m + 1
    end while
    return t * s
end function
printf(1,"    Values of the ?2 probablility distribution function\n")
printf(1," x/k    1         2         3         4         5\n")
for x=0 to 10 do
    printf(1,"%2d  ", x)
    for k=1 to 5 do
        printf(1,"%f  ", pdf(x, k))
    end for
    printf(1,"\n")
end for
printf(1,"\n------Checking cpdf formula works for k = 2-------\n")
printf(1,"    Values of the ?2 cumulative probability distribution function for k = 2\n")
printf(1," x\n")
for x=0 to 10 do
    printf(1,"%2d  %f\n", {x, cpdf(x, 2)})
end for
printf(1,"\n    Values of the ?2 cumulative pdf using simple formula for k = 2\n")
printf(1," x\n")
for x=0 to 10 do
    printf(1,"%2d  %f\n", {x, 1 - exp(-x/2)})
end for

Output:

Same output (so probably not very helpful)

    Values of the ?2 probablility distribution function
 x/k    1         2         3         4         5
 0  0.000000  0.000000  0.000000  0.000000  0.000000
 1  0.241971  0.303265  0.241971  0.151633  0.080657
 2  0.103777  0.183940  0.207554  0.183940  0.138369
 3  0.051393  0.111565  0.154180  0.167348  0.154180
 4  0.026995  0.067668  0.107982  0.135335  0.143976
 5  0.014645  0.041042  0.073225  0.102606  0.122042
 6  0.008109  0.024894  0.048652  0.074681  0.097304
 7  0.004553  0.015099  0.031873  0.052845  0.074371
 8  0.002583  0.009158  0.020667  0.036631  0.055112
 9  0.001477  0.005554  0.013296  0.024995  0.039887
10  0.000850  0.003369  0.008500  0.016845  0.028335

------Checking cpdf formula works for k = 2-------
    Values of the ?2 cumulative probability distribution function for k = 2
 x
 0  0.000000
 1  0.393469
 2  0.632121
 3  0.776870
 4  0.864665
 5  0.917915
 6  0.950213
 7  0.969803
 8  0.981684
 9  0.988891
10  0.993262

    Values of the ?2 cumulative pdf using simple formula for k = 2
 x
 0  0.000000
 1  0.393469
 2  0.632121
 3  0.776870
 4  0.864665
 5  0.917915
 6  0.950213
 7  0.969803
 8  0.981684
 9  0.988891
10  0.993262