montepython_public
montepython_public copied to clipboard
baudren
Covariance matrix has negative eigenvalues
I am doing MCMC with Planck likelihood + additional custom parameters in the power spectrum (using generate_Pk.py). To get a covariance matrix, I did the follows:
(1) Do a short run:
montepython run --conf default.conf -p default.param -c covmat/base.covmat -o res_mcmc_short -N 1000
Here the base.covmat is the default covariance matrix which comes with montepython. Note that base.covmat only has standard parameters in it, not my custom parameters.
(2) After the short run, a covariance matrix is generated with all parameters. Then I plan to use the new covariance for a long run. However, there are eigenvalues of the new covariance matrix. So at the first chain step, the code stops, with error message
File "/home/wangyi/PublicCode/montepython/montepython/mcmc.py", line 116, in get_new_position
rd.gauss(0, 1)*data.jumping_factor
The error message is because in that line of code one calculate the sqrt(eigenvalue_of_covmat), which cannot proceed when the eigenvalue is negative.
Would you see if this is because I am doing something wrong, or it is a bug in the code? Thanks a lot!
===== Attachments =====
(1) My custom spectrum, and generated covariance matrix (in res_mcmc_short folder).
https://www.dropbox.com/s/1l312qzh2kjobs6/vary_all_mcmc.tar.gz?dl=0
(2) Mathematica notebook to show that the covariance matrix has negative eigenvalue.
https://www.dropbox.com/s/h6313jrzg15t9jd/load_covmat.nb?dl=0
This is probably because the "short" run I tried was too short (1000 steps for Planck likelihood). As an attempt, I changed it into 10000 steps, the covmat is then positive definite.
Other comments/suggestions are welcome.
Hello,
yes, this is most likely the explanation. There is no guarantee that a too short run will produce a meaningful covariance matrix, especially with high dimensions. One way of checking if the covariance matrix has a chance to be meaningul is to look at the posterior distribution (especially the triangle plot). If most 2d diagrams are very patchy, you should run longer.
