Statistics/Chi-squared distribution: Difference between revisions
julia example
(Added Phix translation of the code on the talk page) |
(julia example) |
||
Line 77:
<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}}==
| |||