search

hakaru-dev / hakaru

A probabilistic programming language
BSD 3-Clause "New" or "Revised" License
311 stars 30 forks source link

Chi sq tests #142

Closed mkhattab940 closed 6 years ago

mkhattab940 commented 6 years ago

Starting to add test cases for the standardChiSq() distribution as implemented in stdlib.

Here is an initial test to make sure we are on the right track.

Based on the 6th bullet point here: https://en.wikipedia.org/wiki/Chi-squared_distribution#Relation_to_other_distributions

Relevant output from stack test:

### Failure in: 6:RoundTrip:6:0:t_stdChiSq_to_gamma:1
haskell/Tests/TestTools.hs:130
expected:
gamma(9/1, 1/2)
but got:
q30D <~ normal(+0/1, 1/1)
q31B <~ normal(+0/1, 1/1)
q32z <~ normal(+0/1, 1/1)
q33x <~ normal(+0/1, 1/1)
q34v <~ normal(+0/1, 1/1)
q35t <~ normal(+0/1, 1/1)
q36r <~ normal(+0/1, 1/1)
q37p <~ normal(+0/1, 1/1)
q38n <~ normal(+0/1, 1/1)
q39l <~ normal(+0/1, 1/1)
q310j <~ normal(+0/1, 1/1)
q311h <~ normal(+0/1, 1/1)
q312f <~ normal(+0/1, 1/1)
q313d <~ normal(+0/1, 1/1)
q314b <~ normal(+0/1, 1/1)
q3159 <~ normal(+0/1, 1/1)
q3167 <~ normal(+0/1, 1/1)
q3175 <~ normal(+0/1, 1/1)
return real2prob
         (q30D ^ 2 * (+1/4)
          + q31B ^ 2 * (+1/4)
          + q310j ^ 2 * (+1/4)
          + q311h ^ 2 * (+1/4)
          + q312f ^ 2 * (+1/4)
          + q313d ^ 2 * (+1/4)
          + q314b ^ 2 * (+1/4)
          + q3159 ^ 2 * (+1/4)
          + q3167 ^ 2 * (+1/4)
          + q3175 ^ 2 * (+1/4)
          + q32z ^ 2 * (+1/4)
          + q33x ^ 2 * (+1/4)
          + q34v ^ 2 * (+1/4)
          + q35t ^ 2 * (+1/4)
          + q36r ^ 2 * (+1/4)
          + q37p ^ 2 * (+1/4)
          + q38n ^ 2 * (+1/4)
          + q39l ^ 2 * (+1/4))
Cases: 326  Tried: 282  Errors: 3  Failures: 14
Cases: 326  Tried: 283  Errors: 3  Failures: 14
Cases: 326  Tried: 284  Errors: 3  Failures: 14
Cases: 326  Tried: 285  Errors: 3  Failures: 14
JacquesCarette commented 6 years ago

I am not entirely surprised that fails - it is an 18-dimensional integral, after all. Have you tried smaller examples?

JacquesCarette commented 6 years ago

Also, I would prefer if you used rationals instead of floats in the tests. So 9/1 and 1/2.

mkhattab940 commented 6 years ago

OK, I have fixed it to test a base case (v = 1) of the relation. Looks like it still fails:

### Failure in: 6:RoundTrip:6:0:t_stdChiSq_to_gamma:1
haskell/Tests/TestTools.hs:130
expected:
gamma(1/2, 1/1)
but got:
q305 <~ normal(+0/1, 1/1)
return real2prob(q305 ^ 2) * (1/2)
Cases: 326  Tried: 282  Errors: 2  Failures: 15
Cases: 326  Tried: 283  Errors: 2  Failures: 15
Cases: 326  Tried: 284  Errors: 2  Failures: 15
Cases: 326  Tried: 285  Errors: 2  Failures: 15
JacquesCarette commented 6 years ago

Can you derive this 'by hand', aka in Maple, but on the PDFs? Then it would be worth filing a bug report on the main Hakaru repo.

mkhattab940 commented 6 years ago

I have the following derivation done in maple. Do you want me to add the .mpl as a file somewhere?

    |\^/|     Maple 2015 (X86 64 LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2015
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> PDF_chisq := (n, x) -> x^(1/2*n-1)*exp(-1/2*x)/(2^(1/2*n))/GAMMA(1/2*n); 
                                        (1/2 n - 1)
                                       x            exp(-1/2 x)
                PDF_chisq := (n, x) -> ------------------------
                                         (1/2 n)
                                        2        GAMMA(1/2 n)

> PDF_cX := (n, c, x) -> PDF_chisq(n,x/c)/c;
                                          PDF_chisq(n, x/c)
                   PDF_cX := (n, c, x) -> -----------------
                                                  c

> PDF_gamma := (k, t, x) -> PDF_cX(2*k,1/2*t,x);
                PDF_gamma := (k, t, x) -> PDF_cX(2 k, 1/2 t, x)

> simplify(PDF_gamma(k, t, x));
                                    k
                               (x/t)  exp(- x/t)
                               -----------------
                                  x GAMMA(k)

> quit
memory used=2.7MB, alloc=8.3MB, time=0.10
JacquesCarette commented 6 years ago

Yes, please do add the .mpl source file that produces that output in the repository.

It would be even better if you could do it using the PDFs inside the Maple library that supports Hakaru - but this isn't very well documented. I should try to give you an example.

Unfortunately, I am unlikely to get to this for ~10 days (I am away all of next week).

mkhattab940 commented 6 years ago

Ok, I have modified my maple code to better reflect a formal proof and added it in a file.

I have also writted a stdChiSq_to_exponential test which also failed:

### Failure in: 6:RoundTrip:6:1:t_stdChiSq_to_exponential:0
haskell/Tests/TestTools.hs:130
expected:
exponential = fn alpha prob:
              X <~ uniform(nat2real(0), nat2real(1))
              return -prob2real(nat2prob(1))
                      * prob2real(alpha)
                      * log(real2prob(X))
exponential(1/2)
but got:
X3 <~ uniform(+0/1, +1/1)
return log(real2prob(X3)) * (-1/2)
Cases: 328  Tried: 283  Errors: 2  Failures: 16

### Failure in: 6:RoundTrip:6:1:t_stdChiSq_to_exponential:1
haskell/Tests/TestTools.hs:130
expected:
exponential = fn alpha prob:
              X <~ uniform(nat2real(0), nat2real(1))
              return -prob2real(nat2prob(1))
                      * prob2real(alpha)
                      * log(real2prob(X))
exponential(1/2)
but got:
q307 <~ normal(+0/1, 1/1)
q315 <~ normal(+0/1, 1/1)
return q307 ^ 2 + q315 ^ 2

I have also added a maple proof for this test case as well.

JacquesCarette commented 6 years ago

So the first failure isn't really a failure, in that what it returns is (1) correct, and (2) simpler. In some cases, the 'normal form' of both sides of an equivalence will actually be a 3rd expression. So there you'll need to adjust the test.

The second failure is something that does need to be investigated, as it does seem like something which ought to be feasible.

mkhattab940 commented 6 years ago

OK, so I've messed around with it a bit and I can make the first failure go away, but only by dismantling the function and basically just rewriting it the way that it's shown under but got:

Is that how I should do it then?

JacquesCarette commented 6 years ago

So after some more investigation, it looks like this is going to stay a failure. To verify that these are the same requires a change of variables, which the code currently does not automatically do.

Could you move this test to a (new) file, Roundtrip2, where we will expect things to fail? Then I can merge this in.

mkhattab940 commented 6 years ago

OK, done

mkhattab940 commented 6 years ago

OK, I've corrected the rayleigh_to_stdChiSq tests and added a maple proof. Here's the errors this case is giving:

### Failure in: 6:RoundTrip:6:2:t_rayleigh_to_chiSq:0
haskell/Tests/TestTools.hs:130
expected:
chiSq_iid = fn n nat:
            fn mean real:
            fn stdev prob:
            q <~ plate _ of n: normal(mean, stdev)
            return summate i from 0 to size(q):
                   ((q[i] - mean) * prob2real(1/ stdev)) ^ 2
standardChiSq = fn n nat: chiSq_iid(n, nat2real(0), nat2prob(1))
standardChiSq(2)
but got:
q307 <~ normal(+0/1, 1/1)
q315 <~ normal(+0/1, 1/1)
return q307 ^ 2 + q315 ^ 2
Cases: 330  Tried: 285  Errors: 2  Failures: 18

### Failure in: 6:RoundTrip:6:2:t_rayleigh_to_chiSq:1
haskell/Tests/TestTools.hs:130
expected:
chiSq_iid = fn n nat:
            fn mean real:
            fn stdev prob:
            q <~ plate _ of n: normal(mean, stdev)
            return summate i from 0 to size(q):
                   ((q[i] - mean) * prob2real(1/ stdev)) ^ 2
standardChiSq = fn n nat: chiSq_iid(n, nat2real(0), nat2prob(1))
standardChiSq(2)
but got:
X3 <~ uniform(+0/1, +1/1)
return log(real2prob(X3)) * (-2/1)
Cases: 330  Tried: 286  Errors: 2  Failures: 19
Cases: 330  Tried: 287  Errors: 2  Failures: 19
Cases: 330  Tried: 288  Errors: 2  Failures: 19
Cases: 330  Tried: 289  Errors: 2  Failures: 19