View source for 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:

<syntaxhighlight lang="wren">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)
}</syntaxhighlight>

{{out}}
<pre>
    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
</pre>
As there's some heavy math here, a Julia or Python example to check against would be useful. --[[User:PureFox|PureFox]] ([[User talk: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 [[j:Essays/Chi_Squared_CDF|entry]] on this subject which claims to include an example which matches "calculated results against values from the [[WikiPedia:Handbook_of_Mathematical_Functions|Handbook of Mathematical Functions]] by Abramowitz and Stegun, Table 26.8." --[[User:Rdm|Rdm]] ([[User talk: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)

[[User:PureFox|PureFox]] ([[User talk:PureFox|talk]]) 10:48, 2 October 2022 (UTC)
::: This tells us that that notation in the task description matches that notation in the wikipedia entry. --[[User:Rdm|Rdm]] ([[User talk:Rdm|talk]]) 11:11, 2 October 2022 (UTC)

:: I am right now also having trouble with the cdf formula at the level of the incomplete Gamma I am using. I think I may have to calculate it in 2 or 3 different ways depending on x and k. --[[User:Wherrera|Wherrera]] ([[User talk:Wherrera|talk]]) 18:33, 2 October 2022 (UTC)

:: Well, I gave up getting incomplete gamma to work accurately with just two methods, and am just going to call the one in the BigFloat MPFR library, which I think uses six or more methods for incomplete gamma based on the a and x values. --[[User:Wherrera|Wherrera]] ([[User talk:Wherrera|talk]]) 19:43, 2 October 2022 (UTC)

:: I did find that the two formulas I had to calculate the last summation used 2 different meanings for x  which had things overly complicated and incorrect. Fixed. --[[User:Wherrera|Wherrera]] ([[User talk:Wherrera|talk]]) 01:22, 3 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. --[[User:PureFox|PureFox]] ([[User talk:PureFox|talk]]) 09:48, 3 October 2022 (UTC)