Talk:Statistics/Chi-squared distribution
Chi-squared cumulative probability distribution
I'm wondering whether the formula given in the task description for the above is in fact correct?
According to Wikipedia, for the case k = 2, this formula reduces to F(x; 2) = 1 - exp(-x/2). But for values of x from 0 to 10, the only way I can get this simple formula to agree with the one in the task description is if I halve 'x' first before applying the latter. Here's the Wren code I have so far in case I'm missing something or have done something silly:
import "./math" for Math
import "./fmt" for Fmt
class Chi2 {
static pdf(x, k) {
if (x <= 0) return 0 // this should be in task description
return (-x/2).exp * x.pow(k/2-1) / (2.pow(k/2) * Math.gamma(k/2))
}
static cpdf(x, k) {
x = x / 2 // need to do this to agree with Wikipedia formula for k = 2
var t = (-x).exp * x.pow(k/2) * Math.gamma(k/2) / Math.gamma(k)
var s = 0
var m = 0
var tol = 1e-15 // say
while (true) {
var term = x.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 probablility 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------Checking cpdf formula works for k = 2-------")
System.print("\n Values of the χ2 cumulative probability distribution function for k = 2")
System.print(" x")
for (x in 0..10) {
Fmt.print("$2d $f", x, Chi2.cpdf(x, 2))
}
System.print("\n Values of the χ2 cumulative pdf using simple formula for k = 2")
System.print(" x")
for (x in 0..10) {
Fmt.print("$2d $f", x, 1 - (-x/2).exp)
}
- Output:
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
As there's some heavy math here, a Julia or Python example to check against would be useful. --PureFox (talk) 09:46, 2 October 2022 (UTC)
- I haven't looked closely at this one, yet. But it's also worth asking whether the wikipedia entry is using notation in the way we expect
- That said, there's a small j wiki entry on this subject which claims to include an example which matches "calculated results against values from the Handbook of Mathematical Functions by Abramowitz and Stegun, Table 26.8." --Rdm (talk) 09:59, 2 October 2022 (UTC)
Well, the Wikipedia formula looks correct to me. If we substitute k = 2 in their formula for the pdf (which agrees with the task description so there's no notational difference here) we get (since gamma(1) = 1):
- f(x; 2) = exp(-x/2) / 2
and if we integrate that we get:
- F(x; 2) = 1 - exp(-x/2)
PureFox (talk) 10:48, 2 October 2022 (UTC)
Thanks for correcting the cumulative pdf formula. I'm pleased we didn't have to resort to the mpfr_gamma_inc function as it's not one I've wrapped in Wren-GMP. --PureFox (talk) 09:48, 3 October 2022 (UTC)