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}\,e^{-x}\sum _{m=0}^{\infty }{\frac {x^{m}}{\Gamma (s+m+1)}}$$

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

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

Because series formulas may be subject to accumulated errors from rounding in the frequently used region where x and k are under 10 and near one another, you may 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 to 12 decimal places """
function gamma(x)
    p = [ 0.99999999999980993, 676.5203681218851, -1259.1392167224028,
          771.32342877765313, -176.61502916214059, 12.507343278686905,
          -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 ]
    if x < 0.5
        return π / (sinpi(x) * gamma(1.0 - x))
    else
        x -= 1.0
        t = p[1]
        for i in 1:8
            t += p[i+1] / (x + i)
        end
    end
    w = x + 7.5
    return sqrt(2.0 * π) * w^(x+0.5) * exp(-w) * t
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

""" lower incomplete gamma by series formula with gamma """
function gamma_cdf(k, x)
    return x^k * exp(-x) * sum(x^m / gamma(k + m + 1) for m in 0:100)
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

""" 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     cdf for χ2   P value (df=3)\n", "-"^36)
for p in [1, 2, 4, 8, 16, 32]
    cdf = round(cdf_χ2(p, 3), digits=10)
    println(lpad(p, 2), "     ", cdf, "   ",  round(1.0 - cdf, digits=10))
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))

using Plots
x = 0.0:0.01:10
y = [map(p -> χ2(p, k), x) for k in 0:3]

plot(x, y, yaxis=[-0.1, 0.5], labels=[0 1 2 3])

Output:

• 

  Graph of Chi-Squared for k values 0 through 3

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     cdf for χ2   P value (df=3)
------------------------------------
 1     0.1987480431   0.8012519569
 2     0.4275932955   0.5724067045
 4     0.7385358701   0.2614641299
 8     0.9539882943   0.0460117057
16     0.9988660157   0.0011339843
32     0.9999994767   5.233e-7

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

Phix

Translation of: Julia

--
-- demo\rosetta\Chi-squared_distribution.exw
--
without js          -- mpfr_gamma[_inc]()
requires("1.0.2")   -- mpfr_gamma_inc()
include mpfr.e

function gamma(atom x)
-- gamma function for x, using the MPFR library
    mpfr {z, bigx} = mpfr_inits(2,{0,x})
    mpfr_gamma(z, bigx, 0)
    return mpfr_get_d(z)
end function

function chi_squared(atom x, k)
-- Chi-squared function, the probability distribution function (pdf) for chi-squared
    return iff(x > 0 ? power(x,k/2-1) * exp(-x/2) / (power(2,k/2) * gamma(k / 2)) : 0)
end function

function gamma_cdf(atom a, x)
-- gamma CDF from lower incomplete gamma function, using MPFR library to get upper incomplete gamma
    mpfr {z, biga, bigx} = mpfr_inits(3,{0,a,x})
    mpfr_gamma_inc(z, biga, bigx, 0)
    return 1-mpfr_get_d(z)/gamma(a)
end function

function cdf_chi_squared(atom x, k)
-- Cumulative probability function (cdf) for chi-squared
    return iff(x<=0 or k<=0 ? 0.0 : gamma_cdf(k/2, x/2))
end function

printf(1,"       ------------------------------------ Chi-squared ------------------------------------\n")
printf(1," x             k = 1             k = 2             k = 3             k = 4             k = 5\n")
printf(1,repeat('-',93)&"\n")
for x=0 to 10 do
    printf(1,"%2d",x)
    for k=1 to 5 do
        printf(1,"%18.11f%n",{chi_squared(x, k),k=5})
    end for
end for

printf(1,"\nChi_squared x     P value (df=3)\n------------------------------------\n")
for p in {1, 2, 4, 8, 16, 32} do
    printf(1,"      %2d          %.16g\n",{p, 1-cdf_chi_squared(p, 3)})
end for

constant 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 },

fmt = "\n"&"""
For the airport data, diff total is %.15f,
              degrees of freedom is %d,
                      ch-squared is %.15f, 
                         p value is %.15f
"""
integer df = length(airportdata)/2-1
atom dtotal = sum(sq_div(sq_power(sq_sub(airportdata,expected_data),2),expected_data))
printf(1,fmt,{dtotal, df, chi_squared(dtotal,df), cdf_chi_squared(dtotal, df)})

include IupGraph.e

function get_data(Ihandle /*graph*/)
    constant colours = {CD_BLUE,CD_ORANGE,CD_GREEN,CD_RED,CD_PURPLE}
    sequence x = sq_div(tagset(999,0),100),
             xy = {{"NAMES",{"0","1","2","3","4"}}}
    for k=0 to 4 do
        xy = append(xy,{x,apply(true,chi_squared,{x,k}),colours[k+1]})
    end for
    return xy
end function

IupOpen()
Ihandle graph = IupGraph(get_data,`RASTERSIZE=340x180,GRID=NO`)
IupSetAttributes(graph,`XMAX=10,XTICK=2,XMARGIN=10,YMAX=0.5,YTICK=0.1`)
IupShow(IupDialog(graph,`TITLE="Chi-squared distribution",MINSIZE=260x200`))
if platform()!=JS then
    IupMainLoop()
    IupClose()
end if

Output:

       ------------------------------------ Chi-squared ------------------------------------
 x             k = 1             k = 2             k = 3             k = 4             k = 5
---------------------------------------------------------------------------------------------
 0     0.00000000000     0.00000000000     0.00000000000     0.00000000000     0.00000000000
 1     0.24197072452     0.30326532986     0.24197072452     0.15163266493     0.08065690817
 2     0.10377687436     0.18393972059     0.20755374871     0.18393972059     0.13836916581
 3     0.05139344327     0.11156508007     0.15418032980     0.16734762011     0.15418032980
 4     0.02699548326     0.06766764162     0.10798193303     0.13533528324     0.14397591070
 5     0.01464498256     0.04104249931     0.07322491281     0.10260624828     0.12204152135
 6     0.00810869555     0.02489353418     0.04865217333     0.07468060255     0.09730434666
 7     0.00455334292     0.01509869171     0.03187340045     0.05284542099     0.07437126772
 8     0.00258337317     0.00915781944     0.02066698535     0.03663127778     0.05511196094
 9     0.00147728280     0.00555449827     0.01329554524     0.02499524221     0.03988663571
10     0.00085003666     0.00336897350     0.00850036660     0.01684486750     0.02833455534

Chi_squared x     P value (df=3)
------------------------------------
       1          0.8012519569012008
       2          0.5724067044708798
       4          0.2614641299491106
       8          0.0460117056892314
      16          0.0011339842897853
      32          5.233466447984724e-7

For the airport data, diff total is 4.512820512820513,
              degrees of freedom is 3,
                      ch-squared is 0.088753925984435,
                         p value is 0.788850426319306

Python

''' rosettacode.org/wiki/Statistics/Chi-squared_distribution#Python '''


from math import exp, pi, sin, sqrt
from matplotlib.pyplot import plot, legend, ylim


def gamma(x):
    ''' gamma function, accurate to about 12 decimal places '''
    p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
         771.32342877765313, -176.61502916214059, 12.507343278686905,
         -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7]
    if x < 0.5:
        return pi / (sin(pi * x) * gamma(1.0 - x))
    x -= 1.0
    t = p[0]
    for i in range(1, 9):
        t += p[i] / (x + i)

    w = x + 7.5
    return sqrt(2.0 * pi) * w**(x+0.5) * exp(-w) * t


def χ2(x, k):
    ''' Chi-squared function, the probability distribution function (pdf) for chi-squared '''
    return x**(k/2 - 1) * exp(-x/2) / (2**(k/2) * gamma(k / 2)) if x > 0 and k > 0 else 0.0


def gamma_cdf(k, x):
    ''' lower incomplete gamma by series formula with gamma '''
    return x**k * exp(-x) * sum(x**m / gamma(k + m + 1) for m in range(100))


def cdf_χ2(x, k):
    ''' Cumulative probability function (cdf) for chi-squared '''
    return gamma_cdf(k / 2, x / 2) if x > 0 and k > 0 else 0.0


print('x         χ2 k = 1           k = 2           k = 3           k = 4           k = 5')
print('-' * 93)
for x in range(11):
    print(f'{x:2}', end='')
    for k in range(1, 6):
        print(f'{χ2(x, k):16.8}', end='\n' if k % 5 == 0 else '')


print('\nχ2 x     P value (df=3)\n----------------------')
for p in [1, 2, 4, 8, 16, 32]:
    print(f'{p:2}', '    ', 1.0 - cdf_χ2(p, 3))


AIRPORT_DATA = [[77, 23], [88, 12], [79, 21], [81, 19]]

EXPECTED = [[81.25, 18.75],
            [81.25, 18.75],
            [81.25, 18.75],
            [81.25, 18.75]]

DTOTAL = sum((d[pos] - EXPECTED[i][pos])**2 / EXPECTED[i][pos]
             for i, d in enumerate(AIRPORT_DATA) for pos in [0, 1])

print(
    f'\nFor the airport data, diff total is {DTOTAL}, χ2 is {χ2(DTOTAL, 3)}, p value {cdf_χ2(DTOTAL, 3)}')
X = [x * 0.001 for x in range(10000)]
for k in range(5):
    plot(X, [χ2(p, k) for p in X])
legend([0, 1, 2, 3, 4])
ylim(-0.02, 0.5)

Output:

x         χ2 k = 1           k = 2           k = 3           k = 4           k = 5
---------------------------------------------------------------------------------------------
 0             0.0             0.0             0.0             0.0             0.0
 1      0.24197072      0.30326533      0.24197072      0.15163266     0.080656908
 2      0.10377687      0.18393972      0.20755375      0.18393972      0.13836917
 3     0.051393443      0.11156508      0.15418033      0.16734762      0.15418033
 4     0.026995483     0.067667642      0.10798193      0.13533528      0.14397591
 5     0.014644983     0.041042499     0.073224913      0.10260625      0.12204152
 6    0.0081086956     0.024893534     0.048652173     0.074680603     0.097304347
 7    0.0045533429     0.015098692       0.0318734     0.052845421     0.074371268
 8    0.0025833732    0.0091578194     0.020666985     0.036631278     0.055111961
 9    0.0014772828    0.0055544983     0.013295545     0.024995242     0.039886636
10   0.00085003666    0.0033689735    0.0085003666     0.016844867     0.028334555

χ2 x     P value (df=3)
----------------------
 1      0.8012519569012009
 2      0.5724067044708798
 4      0.26146412994911117
 8      0.04601170568923141
16      0.0011339842897852837
32      5.233466447984725e-07

For the airport data, diff total is 4.512820512820513, χ2 is 0.088753925984435, p value 0.7888504263193064

Wren

Library: Wren-math

Library: Wren-fmt

import "./math" for Math
import "./fmt" for Fmt

class Chi2 {
    static pdf(x, k) {
        if (x <= 0) return 0
        return (-x/2).exp * x.pow(k/2-1) / (2.pow(k/2) * Math.gamma(k/2))
    }

    static cpdf(x, k) {
        var t = (-x/2).exp * (x/2).pow(k/2)
        var s = 0
        var m = 0
        var tol = 1e-15 // say
        while (true) {
            var term = (x/2).pow(m) / Math.gamma(k/2 + m + 1)
            s = s + term
            if (term.abs < tol) break
            m = m + 1
        }
        return t * s
    }
}

System.print("    Values of the χ2 probability distribution function")
System.print(" x/k    1         2         3         4         5")
for (x in 0..10) {
    Fmt.write("$2d  ", x)
    for (k in 1..5) {
        Fmt.write("$f  ", Chi2.pdf(x, k))
    }
    System.print()
}

System.print("\n    Values for χ2 with 3 degrees of freedom")
System.print("χ2  cum cpdf  p-value")
for (x in [1, 2, 4, 8, 16, 32]) {
    var cpdf = Chi2.cpdf(x, 3)
    Fmt.print("$2d  $f  $f", x,  cpdf, 1-cpdf)
}

var airport = [[77, 23], [88, 12], [79, 21], [81, 19]]
var expected = [81.25, 18.75]
var dsum = 0
for (i in 0...airport.count) {
    for (j in 0...airport[0].count) {
        dsum = dsum + (airport[i][j] - expected[j]).pow(2) / expected[j]
    }
}
var dof = (airport.count - 1) / (airport[0].count - 1)
System.print("\nFor airport data table: ")
Fmt.print("  diff sum : $f", dsum)
Fmt.print("  d.o.f.   : $d", dof)
Fmt.print("  χ2 value : $f", Chi2.pdf(dsum, 3))
Fmt.print("  p-value  : $f", Chi2.cpdf(dsum, 3))

Output:

    Values of the χ2 probability 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  

    Values for χ2 with 3 degrees of freedom
χ2  cum cpdf  p-value
 1  0.198748  0.801252
 2  0.427593  0.572407
 4  0.738536  0.261464
 8  0.953988  0.046012
16  0.998866  0.001134
32  0.999999  0.000001

For airport data table: 
  diff sum : 4.512821
  d.o.f.   : 3
  χ2 value : 0.088754
  p-value  : 0.788850