BennMacdonald / AGM_RPackage

1 stars 0 forks source link

Check likelihood is calculated correctly #13

Closed BennMacdonald closed 7 years ago

BennMacdonald commented 7 years ago

As you mentioned, since you switched over to a Cholesky decomposition, we must make sure no bug has crept in. I can check this. Once we have the version we wish to upload (after the vignette is written, in case we want to make other changes), I will run the new script for a chosen seed with the previous version I have. I will then check the likelihoods to see if they are (approximately) the same.

FrankD commented 7 years ago

Thanks Ben! I don't think we can test it quite like that, because even with the same random seeds, the two versions of the code will have different sampling paths. This is because every so often, the inversion via solve() in the old version generated a non-positive-definite inverted matrix, resulting in likelihoods that were slightly wrong (or in extreme cases, very wrong, which is why I changed it).

The way I thought to test it is to test the likelihood function in isolation; i.e. generate some test states from a run of the code (either version), and then see if the likelihoods calculated for those states are similar under the old and new version of the function that calculates the likelihood, or if they are different, check if it is due to numerical errors in the solve() inversion.

FrankD commented 7 years ago

I've now checked this, the new likelihood function is correct. However, if you compare it with the old function, you will see non-negligible differences in the log likelihood. When I delved deeper, it became clear that this is simply because solve(X)%*%b is less accurate than backsolve(chol(X), b, transpose=TRUE), i.e. if you front-multiply both by X, the second expression is closer to b by orders of magnitude.

I have now deleted the old likelihood functions to keep the code tidy. I have also added a unit test so that we're alerted if something we do changes the likelihood calculation (for a few random test cases on the LV ODE system).