Statistics/Chi-squared distribution: Difference between revisions
Statistics/Chi-squared distribution (view source)
Revision as of 10:26, 10 February 2024
, 3 months ago→{{header|Wren}}: Changed to Wren S/H
m (gamma) |
m (→{{header|Wren}}: Changed to Wren S/H) |
||
| (27 intermediate revisions by 9 users not shown) | |||
Line 75:
* 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.
;Related Tasks:
* [[Statistics/Basic]]
* [[Statistics/Normal_distribution]]
;See also:
[[https://en.wikipedia.org/wiki/Chi-squared_test Chi Squared Test]]
[[https://www.itl.nist.gov/div898/handbook/eda/section3/eda3666.htm NIST page on the Chi-Square Distribution]]
<br /><br /><br />
=={{header|C++}}==
<syntaxhighlight lang="c++">
#include <iostream>
#include <cmath>
#include <numbers>
#include <iomanip>
#include <array>
// The normalised lower incomplete gamma function.
double gamma_cdf(const double aX, const double aK) {
double result = 0.0;
for ( uint32_t m = 0; m <= 99; ++m ) {
result += pow(aX, m) / tgamma(aK + m + 1);
}
result *= pow(aX, aK) * exp(-aX);
return result;
}
// The cumulative probability function of the Chi-squared distribution.
double cdf(const double aX, const double aK) {
if ( aK > 1'000 && aK < 100 ) {
return 1.0;
}
return ( aX > 0.0 && aK > 0.0 ) ? gamma_cdf(aX / 2, aK / 2) : 0.0;
}
// The probability density function of the Chi-squared distribution.
double pdf(const double aX, const double aK) {
return ( aX > 0.0 ) ? pow(aX, aK / 2 - 1) * exp(-aX / 2) / ( pow(2, aK / 2) * tgamma(aK / 2) ) : 0.0;
}
int main() {
std::cout << " Values of the Chi-squared probability distribution function" << std::endl;
std::cout << " x/k 1 2 3 4 5" << std::endl;
for ( uint32_t x = 0; x <= 10; x++ ) {
std::cout << std::setw(2) << x;
for ( uint32_t k = 1; k <= 5; ++k ) {
std::cout << std::setw(10) << std::fixed << pdf(x, k);
}
std::cout << std::endl;
}
std::cout << "\n Values for the Chi-squared distribution with 3 degrees of freedom" << std::endl;
std::cout << "Chi-squared cumulative pdf p-value" << std::endl;
for ( uint32_t x : { 1, 2, 4, 8, 16, 32 } ) {
const double cumulative_pdf = cdf(x, 3);
std::cout << std::setw(6) << x << std::setw(19) << std::fixed << cumulative_pdf
<< std::setw(14) << ( 1.0 - cumulative_pdf ) << std::endl;
}
const std::array<const std::array<int32_t, 2>, 4> observed =
{ { { 77, 23 }, { 88, 12 }, { 79, 21 }, { 81, 19 } } };
const std::array<const std::array<double, 2>, 4> expected =
{ { { 81.25, 18.75 }, { 81.25, 18.75 }, { 81.25, 18.75 }, { 81.25, 18.75 } } };
double testStatistic = 0.0;
for ( uint64_t i = 0; i < observed.size(); ++i ) {
for ( uint64_t j = 0; j < observed[0].size(); ++j ) {
testStatistic += pow(observed[i][j] - expected[i][j], 2.0) / expected[i][j];
}
}
const uint64_t degreesFreedom = ( observed.size() - 1 ) / ( observed[0].size() - 1 );
std::cout << "\nFor the airport data:" << std::endl;
std::cout << " test statistic : " << std::fixed << testStatistic << std::endl;
std::cout << " degrees of freedom : " << degreesFreedom << std::endl;
std::cout << " Chi-squared : " << std::fixed << pdf(testStatistic, degreesFreedom) << std::endl;
std::cout << " p-value : " << std::fixed << cdf(testStatistic, degreesFreedom) << std::endl;
}
</syntaxhighlight>
{{ out }}
<pre>
Values of the Chi-squared 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 the Chi-squared distribution with 3 degrees of freedom
Chi-squared cumulative pdf 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 the airport data:
test statistic : 4.512821
degrees of freedom : 3
Chi-squared : 0.088754
p-value : 0.788850
</pre>
=={{header|FreeBASIC}}==
{{trans|Wren}}
<syntaxhighlight lang="vb">#define pi 4 * Atn(1)
Function gamma (x As Double) As Double
Dim As Integer k
Dim As Double p(12)
Dim As Double accm = p(1)
If accm = 0 Then
accm = Sqr(2 * pi)
p(1) = accm
Dim As Double k1factrl = 1
For k = 2 To 12
p(k) = Exp(13 - k) * (13 - k)^(k - 1.5) / k1factrl
k1factrl *= -(k - 1)
Next k
End If
For k = 2 To 12
accm += p(k) / (x + k - 1)
Next
accm *= Exp(-(x + 12)) * (x + 12)^(x + .5)
Return accm / x
End Function
Function pdf(x As Double, k As Double) As Double
'probability density function
Return Iif(x <= 0, 0, x^(k / 2 - 1) * Exp(-x / 2) / (2^(k / 2) * gamma(k / 2)))
End Function
Function cpdf(x As Double, k As Double) As Double
'cumulative probability distribution function
Dim As Double t = Exp(-x / 2) * (x / 2)^(k / 2)
Dim As Double s, term
Dim As Uinteger m = 0
Do
term = (x / 2)^m / gamma(k / 2 + m + 1)
s += term
m += 1
Loop Until Abs(term) < 1e-15
Return t * s
End Function
Dim As Integer x, k, i, j
Print " Values of the Chi-squared probability distribution function"
Print " x k = 1 k = 2 k = 3 k = 4 k = 5"
For x = 0 To 10
Print Using "## "; x;
For k = 1 To 5
Print Using "#.############ "; pdf(x, k);
Next
Print
Next
Print !"\n Values for Chi-squared with 3 degrees of freedom"
Print "Chi-squared cum pdf P value"
Dim As Uinteger tt(5) = {1, 2, 4, 8, 16, 32}
For x = 0 To Ubound(tt)
Dim As Double cpdff = cpdf(tt(x), 3)
Print Using " ## #.############ #.############"; tt(x); cpdff; 1-cpdff
Next
Dim As Uinteger airport(3,1) = {{77, 23}, {88, 12}, {79, 21}, {81, 19}}
Dim As Double expected(1) = {81.25, 18.75}
Dim As Double dsum = 0
For i = 0 To Ubound(airport,1)
For j = 0 To Ubound(airport,2)
dsum += (airport(i, j) - expected(j))^2 / expected(j)
Next
Next
Dim As Double dof = Ubound(airport,1) / Ubound(airport,2)
Print Using !"\nFor the airport data, diff total is #.############"; dsum
Print Spc(14); "degrees of freedom is"; dof
Print Spc(21); Using "Chi-squared is #.############"; pdf(dsum, dof)
Print Spc(25); Using "P value is #.############"; cpdf(dsum, dof)
Sleep</syntaxhighlight>
{{out}}
<pre> Values of the Chi-squared probability distribution function
x k = 1 k = 2 k = 3 k = 4 k = 5
0 0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000
1 0.241970724519 0.303265329856 0.241970724519 0.151632664928 0.080656908173
2 0.103776874355 0.183939720586 0.207553748710 0.183939720586 0.138369165807
3 0.051393443268 0.111565080074 0.154180329804 0.167347620111 0.154180329804
4 0.026995483257 0.067667641618 0.107981933026 0.135335283237 0.143975910702
5 0.014644982562 0.041042499312 0.073224912810 0.102606248280 0.122041521349
6 0.008108695555 0.024893534184 0.048652173330 0.074680602552 0.097304346659
7 0.004553342922 0.015098691711 0.031873400451 0.052845420989 0.074371267720
8 0.002583373169 0.009157819444 0.020666985354 0.036631277777 0.055111960944
9 0.001477282804 0.005554498269 0.013295545236 0.024995242211 0.039886635707
10 0.000850036660 0.003368973500 0.008500366603 0.016844867498 0.028334555342
Values for Chi-squared with 3 degrees of freedom
Chi-squared cum pdf P value
1 0.198748043099 0.801251956901
2 0.427593295529 0.572406704471
4 0.738535870051 0.261464129949
8 0.953988294311 0.046011705689
16 0.998866015710 0.001133984290
32 0.999999476654 0.000000523346
For the airport data, diff total is 4.512820512821
degrees of freedom is 3
Chi-squared is 0.088753925984
P value is 0.788850426319</pre>
=={{header|Java}}==
<syntaxhighlight lang="java">
import java.util.List;
public final class StatisticsChiSquaredDistribution {
public static void main(String[] aArgs) {
System.out.println(" Values of the Chi-squared probability distribution function");
System.out.println(" x/k 1 2 3 4 5");
for ( int x = 0; x <= 10; x++ ) {
System.out.print(String.format("%2d", x));
for ( int k = 1; k <= 5; k++ ) {
System.out.print(String.format("%10.6f", pdf(x, k)));
}
System.out.println();
}
System.out.println();
System.out.println(" Values for the Chi-squared distribution with 3 degrees of freedom");
System.out.println("Chi-squared cumulative pdf p-value");
for ( int x : List.of( 1, 2, 4, 8, 16, 32 ) ) {
final double cdf = cdf(x, 3);
System.out.println(String.format("%6d%19.6f%14.6f", x, cdf, ( 1.0 - cdf )));
}
final int[][] observed = { { 77, 23 }, { 88, 12 }, { 79, 21 }, { 81, 19 } };
final double[][] expected = { { 81.25, 18.75 }, { 81.25, 18.75 }, { 81.25, 18.75 }, { 81.25, 18.75 } };
double testStatistic = 0.0;
for ( int i = 0; i < observed.length; i++ ) {
for ( int j = 0; j < observed[0].length; j++ ) {
testStatistic += Math.pow(observed[i][j] - expected[i][j], 2.0) / expected[i][j];
}
}
final int degreesFreedom = ( observed.length - 1 ) / ( observed[0].length - 1 );
System.out.println();
System.out.println("For the airport data:");
System.out.println(" test statistic : " + String.format("%.6f", testStatistic));
System.out.println(" degrees of freedom : " + degreesFreedom);
System.out.println(" Chi-squared : " + String.format("%.6f", pdf(testStatistic, degreesFreedom)));
System.out.println(" p-value : " + String.format("%.6f", cdf(testStatistic, degreesFreedom)));
}
// The gamma function.
private static double gamma(double aX) {
if ( aX < 0.5 ) {
return Math.PI / ( Math.sin(Math.PI * aX) * gamma(1.0 - aX) );
}
final double[] probabilities = new double[] {
0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059,
12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 };
aX -= 1.0;
double t = probabilities[0];
for ( int i = 1; i < 9; i++ ) {
t += probabilities[i] / ( aX + i );
}
final double w = aX + 7.5;
return Math.sqrt(2.0 * Math.PI) * Math.pow(w, aX + 0.5) * Math.exp(-w) * t;
}
// The probability density function of the Chi-squared distribution.
private static double pdf(double aX, double aK) {
return ( aX > 0.0 ) ?
Math.pow(aX, aK / 2 - 1) * Math.exp(-aX / 2) / ( Math.pow(2, aK / 2) * gamma(aK / 2) ) : 0.0;
}
// The cumulative probability function of the Chi-squared distribution.
private static double cdf(double aX, double aK) {
if ( aX > 1_000 && aK < 100 ) {
return 1.0;
}
return ( aX > 0.0 && aK > 0.0 ) ? gammaCDF(aX / 2, aK / 2) : 0.0;
}
// The normalised lower incomplete gamma function.
private static double gammaCDF(double aX, double aK) {
double result = 0.0;
for ( int m = 0; m <= 99; m++ ) {
result += Math.pow(aX, m) / gamma(aK + m + 1);
}
result *= Math.pow(aX, aK) * Math.exp(-aX);
return result;
}
}
</syntaxhighlight>
{{ out }}
<pre>
Values of the Chi-squared 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 the Chi-squared distribution with 3 degrees of freedom
Chi-squared cumulative pdf 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 the airport data:
test statistic : 4.512821
degrees of freedom : 3
Chi-squared : 0.088754
p-value : 0.788850
</pre>
=={{header|jq}}==
{{works with|jq}}
'''Also works with gojq, the Go implementation of jq.'''
The implementation of Chi_cdf here uses the recusive relation for the
gamma function for efficiency and to simplify the avoidance of
numerical overflow issues.
'''Formatting'''
<syntaxhighlight lang=jq>
def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;
# For formatting numbers
# ... but leave non-numbers, 0, and numbers using E alone
def round($dec):
def rpad($len): tostring | ($len - length) as $l | . + ("0" * $l);
if type == "string" then .
elif . == 0 then "0"
else pow(10;$dec) as $m
| . * $m | round / $m
| tostring
| (capture("(?<n>^[^.]*[.])(?<f>[0-9]*$)")
| .n + (.f|rpad($dec)))
// if test("^[0-9]+$") then . + "." + ("" | rpad($dec)) else null end
// .
end;
</syntaxhighlight>
'''Chi-squared pdf and cdf'''
<syntaxhighlight lang=jq>
def Chi2_pdf($x; $k):
if $x <= 0 then 0
else # ((-$x/2)|exp) * pow($x; $k/2 -1) / (pow(2;$k/2) * (($k/2)|gamma))
((-$x/2) + ($k/2 -1) * ($x|log) - ( ($k/2)*(2|log) + (($k/2)|lgamma))) | exp
end;
# $k is the degrees of freedom
# Use recursive relation for gamma: G(x+1) = x * G(x)
def Chi2_cdf($x; $k):
if $x == 0 then 0
elif $x > (1e3 * $k) then 1
else 1e-15 as $tol # for example
| { s: 0, m: 0, term: (1 / ((($k/2)+1)|gamma)) }
| until (.term|length < $tol;
.s += .term
| .m += 1
| .term *= (($x/2) / (($k/2) + .m )) )
| .s * ( ((-$x/2) + ($k/2)*(($x/2)|log)) | exp)
end ;
</syntaxhighlight>
'''The Tasks'''
<syntaxhighlight lang=jq>
def tables:
def r: round(6) | lpad(10);
" Values of the χ2 probability distribution function",
" x/k 1 2 3 4 5",
(range(0;11) as $x
| "\($x|lpad(2)) \(reduce range(1; 6) as $k (""; . + (Chi2_pdf($x; $k) | r) ))" ),
"\n Values for χ2 with 3 degrees of freedom",
"χ2 cum pdf p-value",
( (1, 2, 4, 8, 16, 32) as $x
| Chi2_cdf($x; 3) as $cdf
| "\($x|lpad(2)) \($cdf|r)\(1-$cdf|r)"
);
def airport: [[77, 23], [88, 12], [79, 21], [81, 19]];
def expected: [81.25, 18.75];
def airport($airport; $expected):
def dof: ($airport|length - 1) / ($airport[0]|length - 1);
reduce range(0; $airport|length) as $i (0;
reduce range(0; $airport[0]|length) as $j (.;
. + (($airport[$i][$j] - $expected[$j]) | .*.) / $expected[$j] ) )
| ("\nFor airport data table:",
" diff sum : \(.)",
" d.o.f. : \(dof)",
" χ2 value : \(Chi2_pdf(.; dof))",
" p-value : \(Chi2_cdf(.; dof)|round(4))" ) ;
tables,
airport(airport; expected)
</syntaxhighlight>
{{output}}
<pre>
Values of the χ2 probability distribution function
x/k 1 2 3 4 5
0 0 0 0 0 0
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 pdf 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 1e-06
For airport data table:
diff sum : 4.512820512820513
d.o.f. : 3
χ2 value : 0.088753925984435
p-value : 0.7889
</pre>
=={{header|Julia}}==
Line 82 ⟶ 538:
""" gamma function to 12 decimal places """
function gamma(
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))
Line 114 ⟶ 569:
function cdf_χ2(x, k)
return x <= 0 || k <= 0 ? 0.0 : gamma_cdf(k / 2, x / 2)
end
Line 131 ⟶ 581:
end
println("\nχ2 x cdf for χ2 P value (df=3)\n", "-"^36)
for p in [1, 2, 4, 8, 16, 32]
println(lpad(p, 2), " ", cdf, " ", round(1.0 - cdf, digits=10))
end
Line 175 ⟶ 625:
10 0.00085003666 0.00336897349 0.00850036660 0.01684486749 0.02833455534
χ2 x cdf for χ2 P value (df=3)
------------------------------------
1 0.
2 0.
4 0.
8 0.
16 0.
32 0.9999994767 5.
For the airport data, diff total is 4.512820512820512, χ2 is 0.08875392598443503, p value 0.7888504263193064
</pre>
=={{header|Nim}}==
{{libheader|gnuplotlib.nim}}
Adapted from Python solution.
<syntaxhighlight lang="Nim">import std/[math, strformat, sugar]
func χ²(x: float; k: int): float =
## χ² function, the probability distribution function (pdf) for χ².
if x > 0: x.pow(k / 2 - 1) * exp(-x/2) / (2.pow(k / 2) * gamma(k / 2)) else: 0
func gammaCdf(k, x: float): float =
## Lower incomplete gamma by series formula with gamma.
for m in 0..99:
result += x^m / gamma(k + m.toFloat + 1)
result *= pow(x, k) * exp(-x)
func cdfχ²(x: float; k: int): float =
## Cumulative probability function (cdf) for χ².
if x > 0 and k > 0: gammaCdf(k / 2, x / 2) else: 0
echo " Values of the χ2 probability distribution function"
echo " x/k 1 2 3 4 5"
for x in 0..10:
stdout.write &"{x:2} "
for k in 1..5:
stdout.write &"{χ²(x.toFloat, k):10.6f}"
echo()
echo()
echo " Values for χ2 with 3 degrees of freedom"
echo "χ² cum pdf p-value"
for x in [1, 2, 4, 8, 16, 32]:
let cdf = cdfχ²(x.toFloat, 3)
echo &"{x:2} {cdf:10.6f} {1 - cdf:10.6f}"
const
AirportData = [[float 77, 23], [float 88, 12], [float 79, 21], [float 81, 19]]
Expected = [[81.25, 18.75], [81.25, 18.75], [81.25, 18.75], [81.25, 18.75]]
var dtotal = 0.0
for pos in [0, 1]:
for i, d in AirportData:
dtotal += (d[pos] - Expected[i][pos])^2 / Expected[i][pos]
echo()
echo &"For the airport data, diff total is {dtotal:.6f}, " &
&"χ² is {χ²(dtotal, 3):.6f}, p value is {cdfχ²(dtotal, 3):.6f}."
### Stretch task ###
import gnuplot
let x = collect(for n in 0..9999: n / 1000)
withGnuplot:
cmd "set multiplot"
cmd "set yrange [-0.02:0.5]"
for k in 0..3:
let y = collect(for p in x: χ²(p, k))
plot(x, y, &"k = {k}", "with lines")
</syntaxhighlight>
{{out}}
<pre> 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
χ² cum pdf 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 the airport data, diff total is 4.512821, χ² is 0.088754, p value is 0.788850.
</pre>
[[File:Chi-squared-nim.png]]
=={{header|Perl}}==
{{trans|Raku}}
<syntaxhighlight lang="perl" line>use v5.36;
use Math::MPFR;
use List::Util 'sum';
sub Gamma ($z) {
my $g = Math::MPFR->new();
Math::MPFR::Rmpfr_gamma($g, Math::MPFR->new($z), 0);
$g;
}
sub chi2($x,$k) { $x>0 && $k>0 ? ($x**($k/2 - 1) * exp(-$x/2)/(2**($k/2)*Gamma($k / 2))) : 0 }
sub gamma_cdf($k,$x) { $x**$k * exp(-$x) * sum map { $x** $_ / Gamma($k+$_+1) } 0..100 }
sub cdf_chi2($x,$k) { ($x <= 0 or $k <= 0) ? 0.0 : gamma_cdf($k / 2, $x / 2) }
print 'x χ² ';
print "k = $_" . ' 'x13 for 1..5;
say "\n" . '-' x (my $width = 93);
for my $x (0..10) {
printf '%2d', $x;
printf ' %.' . (int(($width-2)/5)-4) . 'f', chi2($x, $_) for 1..5;
say '';
}
say "\nχ² x cdf for χ² P value (df=3)\n" . '-' x 36;
for my $p (map { 2**$_ } 0..5) {
my $cdf = cdf_chi2($p, 3);
printf "%2d %-.10f %-.10f\n", $p, $cdf, 1-$cdf;
}
my @airport = <77 23 88 12 79 21 81 19>;
my @expected = split ' ', '81.25 18.75 ' x 4;
my $dtotal;
$dtotal += ($airport[$_] - $expected[$_])**2 / $expected[$_] for 0..$#airport;
printf "\nFor the airport data, diff total is %.5f, χ² is %.5f, p value %.5f\n", $dtotal, chi2($dtotal, 3), cdf_chi2($dtotal, 3);</syntaxhighlight>
{{out}}
<pre>x χ² k = 1 k = 2 k = 3 k = 4 k = 5
---------------------------------------------------------------------------------------------
0 0.00000000000000 0.00000000000000 0.00000000000000 0.00000000000000 0.00000000000000
1 0.24197072451914 0.30326532985632 0.24197072451914 0.15163266492816 0.08065690817305
2 0.10377687435515 0.18393972058572 0.20755374871030 0.18393972058572 0.13836916580686
3 0.05139344326792 0.11156508007421 0.15418032980377 0.16734762011132 0.15418032980377
4 0.02699548325659 0.06766764161831 0.10798193302638 0.13533528323661 0.14397591070183
5 0.01464498256193 0.04104249931195 0.07322491280963 0.10260624827987 0.12204152134939
6 0.00810869555494 0.02489353418393 0.04865217332964 0.07468060255180 0.09730434665928
7 0.00455334292164 0.01509869171116 0.03187340045148 0.05284542098906 0.07437126772012
8 0.00258337316926 0.00915781944437 0.02066698535409 0.03663127777747 0.05511196094425
9 0.00147728280398 0.00555449826912 0.01329554523581 0.02499524221105 0.03988663570744
10 0.00085003666025 0.00336897349954 0.00850036660252 0.01684486749771 0.02833455534173
χ² x cdf for χ² 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 0.0000005233
For the airport data, diff total is 4.51282, χ² is 0.08875, p value 0.78885</pre>
=={{header|Phix}}==
{{trans|Julia}}
{{libheader|Phix/online}}
You can run this online [http://phix.x10.mx/p2js/chisqd.htm here].
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #000080;font-style:italic;">--
-- demo\rosetta\Chi-squared_distribution.exw
--</span>
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">gamma</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">0.99999999999980993</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">676.5203681218851</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">1259.1392167224028</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">771.32342877765313</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">176.61502916214059</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">12.507343278686905</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.13857109526572012</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">9.9843695780195716e-6</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">1.5056327351493116e-7</span> <span style="color: #0000FF;">}</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">z</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0.5</span> <span style="color: #008080;">then</span>
<span style="color: #008080;">return</span> <span style="color: #004600;">PI</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">(</span><span style="color: #7060A8;">sin</span><span style="color: #0000FF;">(</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">*</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">gamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">z</span><span style="color: #0000FF;">))</span>
Line 212 ⟶ 819:
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">
<span style="color: #000080;font-style:italic;">-- Chi-squared function, the probability distribution function (pdf) for chi-squared</span>
<span style="color: #008080;">return</span> <span style="color: #7060A8;">
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">
<span style="color: #004080;">atom</span> <span style="color: #000000;">tot</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">=</span><span style="color: #000000;">
<span style="color: #000000;">tot</span> <span style="color: #0000FF;">+=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">gamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">
<span style="color: #008080;">return</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">tot</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #
<span style="color: #000080;font-style:italic;">-- Cumulative probability function (cdf) for chi-squared</span>
<span style="color: #008080;">return</span> <span style="color: #7060A8;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">0</span> <span style="color: #008080;">or</span> <span style="color: #000000;">k</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">0</span> <span style="color: #0000FF;">?</span> <span style="color: #000000;">0.0</span> <span style="color: #0000FF;">:</span> <span style="color: #000000;">gamma_cdf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">k</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" ------------------------------------ Chi-squared ------------------------------------\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" x k = 1 k = 2 k = 3 k = 4 k = 5\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #008000;">'-'</span><span style="color: #0000FF;">,</span><span style="color: #000000;">92</span><span style="color: #0000FF;">)&</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">to</span> <span style="color: #000000;">10</span> <span style="color: #008080;">do</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%2d
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">5</span> <span style="color: #008080;">do</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\nChi_squared x
<span style="color: #
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">airportdata</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">77</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">88</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">12</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">79</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">21</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">81</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">19</span> <span style="color: #0000FF;">},</span>
<span style="color: #000000;">expected_data</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">81.25</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.75</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">81.25</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.75</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">81.25</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.75</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">81.25</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">18.75</span> <span style="color: #0000FF;">},</span>
<span style="color: #000000;">fmt</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">"\n"</span><span style="color: #0000FF;">&</span><span style="color: #008000;">"""
For the airport data, diff total is %.15f,
degrees of freedom is %d,
ch-squared is %.15f,
p value is %.15f
"""</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">df</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">airportdata</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">dtotal</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_power</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">airportdata</span><span style="color: #0000FF;">,</span><span style="color: #000000;">expected_data</span><span style="color: #0000FF;">),</span><span style="color: #000000;">2</span><span style="color: #0000FF;">),</span><span style="color: #000000;">expected_data</span><span style="color: #0000FF;">))</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">dtotal</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">df</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">chi_squared</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dtotal</span><span style="color: #0000FF;">,</span><span style="color: #000000;">df</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">cdf_chi_squared</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dtotal</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">df</span><span style="color: #0000FF;">)})</span>
<span style="color: #008080;">include</span> <span style="color: #7060A8;">IupGraph</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">get_data</span><span style="color: #0000FF;">(</span><span style="color: #004080;">Ihandle</span> <span style="color: #000080;font-style:italic;">/*graph*/</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">colours</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #004600;">CD_BLUE</span><span style="color: #0000FF;">,</span><span style="color: #004600;">CD_ORANGE</span><span style="color: #0000FF;">,</span><span style="color: #004600;">CD_GREEN</span><span style="color: #0000FF;">,</span><span style="color: #004600;">CD_RED</span><span style="color: #0000FF;">,</span><span style="color: #004600;">CD_PURPLE</span><span style="color: #0000FF;">}</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">tagset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">999</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">),</span><span style="color: #000000;">100</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">xy</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span><span style="color: #008000;">"NAMES"</span><span style="color: #0000FF;">,{</span><span style="color: #008000;">"0"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"1"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"2"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"3"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"4"</span><span style="color: #0000FF;">}}}</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">to</span> <span style="color: #000000;">4</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">xy</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xy</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">apply</span><span style="color: #0000FF;">(</span><span style="color: #004600;">true</span><span style="color: #0000FF;">,</span><span style="color: #000000;">chi_squared</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">}),</span><span style="color: #000000;">colours</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">xy</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #7060A8;">IupOpen</span><span style="color: #0000FF;">()</span>
<span style="color: #004080;">Ihandle</span> <span style="color: #000000;">graph</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">IupGraph</span><span style="color: #0000FF;">(</span><span style="color: #000000;">get_data</span><span style="color: #0000FF;">,</span><span style="color: #008000;">`RASTERSIZE=340x180,GRID=NO`</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">IupSetAttributes</span><span style="color: #0000FF;">(</span><span style="color: #000000;">graph</span><span style="color: #0000FF;">,</span><span style="color: #008000;">`XMAX=10,XTICK=2,XMARGIN=10,YMAX=0.5,YTICK=0.1`</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">IupShow</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">IupDialog</span><span style="color: #0000FF;">(</span><span style="color: #000000;">graph</span><span style="color: #0000FF;">,</span><span style="color: #008000;">`TITLE="Chi-squared distribution",MINSIZE=260x200`</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">platform</span><span style="color: #0000FF;">()!=</span><span style="color: #004600;">JS</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">IupMainLoop</span><span style="color: #0000FF;">()</span>
<span style="color: #7060A8;">IupClose</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
------------------------------------ Chi-squared ------------------------------------
x
--------------------------------------------------------------------------------------------
10 0.00085003666 0.00336897350 0.00850036660 0.01684486750 0.02833455534
Chi_squared x P value (df=3)
------------------------------------
1 0.8012519569012007
2 0.5724067044708797
4 0.2614641299491101
8 0.0460117056892316
16 0.0011339842897863
32 5.233466446874501e-7
For the airport data, diff total is 4.512820512820513,
degrees of freedom is 3,
ch-squared is 0.088753925984435,
p value is 0.788850426319307
</pre>
[[File:Phixplotchisquared.png|thumb|center]]
=={{header|Python}}==
Line 303 ⟶ 931:
from math import
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
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
Line 377 ⟶ 1,024:
For the airport data, diff total is 4.512820512820513, χ2 is 0.088753925984435, p value 0.7888504263193064
</pre>
=={{header|Raku}}==
{{trans|Julia}}
<syntaxhighlight lang="raku" line># 20221101 Raku programming solution
use Graphics::PLplot;
sub Γ(\z) { # https://rosettacode.org/wiki/Gamma_function#Raku
constant g = 9;
z < .5 ?? π/ sin(π * z) / Γ(1 - z)
!! sqrt(2*π) * (z + g - 1/2)**(z - 1/2) * exp(-(z + g - 1/2)) *
[+] < 1.000000000000000174663 5716.400188274341379136
-14815.30426768413909044 14291.49277657478554025
-6348.160217641458813289 1301.608286058321874105
-108.1767053514369634679 2.605696505611755827729
-0.7423452510201416151527e-2 0.5384136432509564062961e-7
-0.4023533141268236372067e-8 > Z* 1, |map 1/(z + *), 0..*
}
sub χ2(\x,\k) {x>0 && k>0 ?? (x**(k/2 - 1)*exp(-x/2)/(2**(k/2)*Γ(k / 2))) !! 0}
sub Γ_cdf(\k,\x) { x**k * exp(-x) * sum( ^101 .map: { x** $_ / Γ(k+$_+1) } ) }
sub cdf_χ2(\x,\k) { (x <= 0 or k <= 0) ?? 0.0 !! Γ_cdf(k / 2, x / 2) }
say ' 𝒙 χ² ', [~] (1..5)».&{ "𝒌 = $_" ~ ' ' x 13 };
say '-' x my \width = 93;
for 0..10 -> \x {
say x.fmt('%2d'), [~] (1…5)».&{χ2(x, $_).fmt: " %-.{((width-2) div 5)-4}f"}
}
say "\nχ² 𝒙 cdf for χ² P value (df=3)\n", '-' x 36;
for 2 «**« ^6 -> \p {
my $cdf = cdf_χ2(p, 3).fmt: '%-.10f';
say p.fmt('%2d'), " $cdf ", (1-$cdf).fmt: '%-.10f'
}
my \airport = [ <77 23>, <88 12>, <79 21>, <81 19> ];
my \expected = [ <81.25 18.75> xx 4 ];
my \dtotal = ( (airport »-« expected)»² »/» expected )».List.flat.sum;
say "\nFor the airport data, diff total is ",dtotal,", χ² is ", χ2(dtotal, 3), ", p value ", cdf_χ2(dtotal, 3);
given Graphics::PLplot.new( device => 'png', file-name => 'output.png' ) {
.begin;
.pen-width: 3 ;
.environment: x-range => [-1.0, 10.0], y-range => [-0.1, 0.5], just => 0 ;
.label: x-axis => '', y-axis => '', title => 'Chi-squared distribution' ;
for 0..3 -> \𝒌 {
.color-index0: 1+2*𝒌;
.line: (0, .1 … 10).map: -> \𝒙 { ( 𝒙, χ2( 𝒙, 𝒌 ) )».Num };
.text-viewport: side=>'t', disp=>-𝒌-2, pos=>.5, just=>.5, text=>'k = '~𝒌
} # plplot.sourceforge.net/docbook-manual/plplot-html-5.15.0/plmtex.html
.end
}</syntaxhighlight>
{{out}}
<pre>
𝒙 χ² 𝒌 = 1 𝒌 = 2 𝒌 = 3 𝒌 = 4 𝒌 = 5
---------------------------------------------------------------------------------------------
0 0.00000000000000 0.00000000000000 0.00000000000000 0.00000000000000 0.00000000000000
1 0.24197072451914 0.30326532985632 0.24197072451914 0.15163266492816 0.08065690817305
2 0.10377687435515 0.18393972058572 0.20755374871030 0.18393972058572 0.13836916580686
3 0.05139344326792 0.11156508007421 0.15418032980377 0.16734762011132 0.15418032980377
4 0.02699548325659 0.06766764161831 0.10798193302638 0.13533528323661 0.14397591070183
5 0.01464498256193 0.04104249931195 0.07322491280963 0.10260624827987 0.12204152134939
6 0.00810869555494 0.02489353418393 0.04865217332964 0.07468060255180 0.09730434665928
7 0.00455334292164 0.01509869171116 0.03187340045148 0.05284542098906 0.07437126772012
8 0.00258337316926 0.00915781944437 0.02066698535409 0.03663127777747 0.05511196094425
9 0.00147728280398 0.00555449826912 0.01329554523581 0.02499524221105 0.03988663570744
10 0.00085003666025 0.00336897349954 0.00850036660252 0.01684486749771 0.02833455534173
χ² 𝒙 cdf for χ² 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 0.0000005233
For the airport data, diff total is 4.512821, χ² is 0.08875392598443506, p value 0.7888504263193072</pre>
[[File:RakuPlotChiSquared.png|thumb|center]]
=={{header|Wren}}==
{{libheader|DOME}}
{{libheader|Wren-math}}
{{libheader|Wren-iterate}}
{{libheader|Wren-fmt}}
{{libheader|Wren-plot}}
<syntaxhighlight lang="wren">import "dome" for Window
import "graphics" for Canvas, Color
import "./math2" for Math
import "./iterate" for Stepped
import "./fmt" for Fmt
import "./plot" for Axes
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 pdf 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, dof))
Fmt.print(" p-value : $f", Chi2.cpdf(dsum, dof))
// generate points for plot
var Pts = List.filled(5, null)
for (k in 0..4) {
Pts[k] = []
var x = 0
while (x < 10) {
Pts[k].add([x, 10 * Chi2.pdf(x, k)])
x = x + 0.01
}
}
class Main {
construct new() {
Window.title = "Chi-squared distribution for k in [0, 4]"
Canvas.resize(1000, 600)
Window.resize(1000, 600)
Canvas.cls(Color.white)
var axes = Axes.new(100, 500, 800, 400, -0.5..10, -0.5..5)
axes.draw(Color.black, 2)
var xMarks = 0..10
var yMarks = 0..5
axes.mark(xMarks, yMarks, Color.black, 2)
var xMarks2 = Stepped.new(0..10, 2)
var yMarks2 = 0..5
axes.label(xMarks2, yMarks2, Color.black, 2, Color.black, 1, 10)
var colors = [Color.blue, Color.yellow, Color.green, Color.red, Color.purple]
for (k in 0..4) {
axes.lineGraph(Pts[k], colors[k], 2)
}
axes.rect(8, 5, 120, 110, Color.black)
for (k in 0..4) {
var y = 4.75 - k * 0.25
axes.line(8.2, y, 9, y, colors[k], 2)
y = 385 - k * 18
axes.print(750, y, k.toString, Color.black)
}
}
init() {}
update() {}
draw(alpha) {}
}
var Game = Main.new()</syntaxhighlight>
{{out}}
Terminal output:
<pre>
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 pdf 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
</pre>
[[File:Wren_chisquared.png|thumb|center]]
| |||