Statistics/Chi-squared distribution: Difference between revisions
Content added Content deleted
(Added Phix translation of the code on the talk page) |
(julia example) |
||
Line 77: | Line 77: | ||
<br /><br /><br /> |
<br /><br /><br /> |
||
=={{header|Julia}}== |
|||
<syntaxhighlight lang="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 |
|||
""" lower incomplete gamma function, by using MPFR ligrary to get upper incomplete gamma """ |
|||
function γ(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) |
|||
q = z/gamma(a) |
|||
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 : γ(k / 2, x / 2) # first(gamma_inc(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, χ2 is ", χ2(dtotal, 3), ", p value ", cdf_χ2(dtotal, 3)) |
|||
</syntaxhighlight>{{out}} |
|||
<pre> |
|||
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, χ2 is 0.08875392598443503, p value 0.7888504263193064 |
|||
</pre> |
|||
=={{header|Phix}}== |
=={{header|Phix}}== |