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.

Julia

""" Rosetta Code task rosettacode.org/wiki/Statistics/Chi-squared_distribution """


""" gamma function for x, using the MPFR library  """
function gamma(x)
    isnan(x) && return x
    z, bigx = BigFloat(), BigFloat(x)
    ccall((:mpfr_gamma, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, 0)
    isnan(z) && throw(DomainError(x, "NaN for gamma in BigFloat library"))
    return Float64(z)
end

""" Chi-squared function, the probability distribution function (pdf) for chi-squared """
function χ2(x, k)
    return x > 0 ? x^(k/2 - 1) * exp(-x/2) / (2^(k/2) * gamma(k / 2)) : 0
end

""" gamma CDF from lower incomplete gamma function, using MPFR library to get upper incomplete gamma """
function gamma_cdf(a, x)
    z, biga, bigx = BigFloat(), BigFloat(a), BigFloat(x)
    ccall((:mpfr_gamma_inc, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32), z, biga, bigx, 0)
    return  Float64(1 - z / gamma(a))
end

""" Cumulative probability function (cdf) for chi-squared """
function cdf_χ2(x, k)
    return x <= 0 || k <= 0 ? 0.0 : gamma_cdf(k / 2, x / 2)  
end

println("x           χ2 k = 1             k = 2             k = 3             k = 4             k = 5")
println("-"^93)
for x in 0:10
      print(lpad(x, 2))
      for k in 1:5
        s = string(χ2(x, k))
        print(lpad(s[1:min(end, 13)], 18), k % 5 == 0 ? "\n" : "")
      end
end


println("\nχ2 x     P value (df=3)\n----------------------")
for p in [1, 2, 4, 8, 16, 32]
    println(lpad(p, 2), "     ", 1.0 - cdf_χ2(p, 3))
end

airportdata = [ 77 23 ;
                88 12;
                79 21;
                81 19 ]

expected_data = [ 81.25 18.75 ;
                  81.25 18.75 ;
                  81.25 18.75 ;
                  81.25 18.75 ; ]

dtotal = sum((airportdata[i] - expected_data[i])^2/ expected_data[i] for i in 1:length(airportdata))

println("\nFor the airport data, diff total is $dtotal, χ2 is ", χ2(dtotal, 3), ", p value ", cdf_χ2(dtotal, 3))

Output:

x           χ2 k = 1             k = 2             k = 3             k = 4             k = 5
---------------------------------------------------------------------------------------------
 0                 0                 0                 0                 0                 0
 1     0.24197072451     0.30326532985     0.24197072451     0.15163266492     0.08065690817
 2     0.10377687435     0.18393972058     0.20755374871     0.18393972058     0.13836916580
 3     0.05139344326     0.11156508007     0.15418032980     0.16734762011     0.15418032980
 4     0.02699548325     0.06766764161     0.10798193302     0.13533528323     0.14397591070
 5     0.01464498256     0.04104249931     0.07322491280     0.10260624827     0.12204152134
 6     0.00810869555     0.02489353418     0.04865217332     0.07468060255     0.09730434665
 7     0.00455334292     0.01509869171     0.03187340045     0.05284542098     0.07437126772
 8     0.00258337316     0.00915781944     0.02066698535     0.03663127777     0.05511196094
 9     0.00147728280     0.00555449826     0.01329554523     0.02499524221     0.03988663570
10     0.00085003666     0.00336897349     0.00850036660     0.01684486749     0.02833455534

χ2 x     P value (df=3)
----------------------
 1     0.8012519569012008
 2     0.5724067044708798
 4     0.2614641299491106
 8     0.04601170568923141
16     0.0011339842897852837
32     5.233466447984725e-7

For the airport data, diff total is 4.512820512820512, χ2 is 0.08875392598443503, p value 0.7888504263193064

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