Using the created covariance matrix

[Kshitizpandey123]

I ran 50000 steps with the Planck-2018 TT data firstly without a covariance matrix and got an abysmal acceptance rate of ~0.01. During this run, it often gave this message ("failed to calculate covariance matrix"), but using the "info" command I anyways generated the matrix.

Now, I am trying to rerun the chains using this covariance matrix as an input (as prescribed in the documentation). But now I am getting this error message at the end: Traceback (most recent call last): File "/home/kshitiz/montepython_public/montepython/MontePython.py", line 42, in sys.exit(run()) File "/home/kshitiz/montepython_public/montepython/run.py", line 45, in run sampler.run(cosmo, data, command_line) File "/home/kshitiz/montepython_public/montepython/sampler.py", line 46, in run mcmc.chain(cosmo, data, command_line) File "/home/kshitiz/montepython_public/montepython/mcmc.py", line 353, in chain Cholesky = la.cholesky(C).T File "/home/kshitiz/.local/lib/python3.10/site-packages/scipy/linalg/_decomp_cholesky.py", line 89, in cholesky c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=True, File "/home/kshitiz/.local/lib/python3.10/site-packages/scipy/linalg/_decomp_cholesky.py", line 37, in _cholesky raise LinAlgError("%d-th leading minor of the array is not positive " numpy.linalg.LinAlgError: 18-th leading minor of the array is not positive definite

It is uploading the matrix correctly, but looks like there is some problem with the matrix itself(?). Pls help.

Regards

[brinckmann]

Hi,

When it's difficult to accumulate points it can be tricky to get that first covmat. I would try to start in a new directory with the example one that's provided for TTonly+lensing and see if that can help you get started, even if a better covmat might need to be computed later (done automatically by the update algorithm): covmat/base2018TTonly_lensing.covmat I'd also recommend to run with the flag --superupdate 20 because it will automatically adjust the jumping factor (starting guess can be provided with the -f float flag), so it will be automatically rescaling the magnitude of the jumps of the MCMC. This means if your starting guesses are poor (which is likely what happened in your first attempt), the code has a better chance to recover than with the covmat update algorithm alone (the latter is activated by default).

If you aren't running with lensing and/or the covmat above doesn't help, I suggest going to the Planck results tables and finding 1-sigma constraints for the cosmological parameters so you can adjust them to your likelihood combination by passing them as the starting guess in the param file (note passing a covmat overrides this number). That's the second to last number (third to last entry), e.g. for tau_reio here it's the 0.008 number: data.parameters['tau_reio'] = [ 0.0543, 0.004, None, 0.008, 1, 'cosmo']

Best, Thejs

[Kshitizpandey123]

Thank you, this helped. I used the provided best-fit values and the covmat/base2018TTonly_lensing.covmat with superupdate. The acceptance rate is now around 0.25 and the plots look nicer too with more data points.

I am trying to use --method Fisher, but it gives these messages during the run: Fisher matrix computation failed - inverse not positive definite

then finally it ended with following messages: numpy.linalg.LinAlgError: 4-th leading minor of the array is not positive definite Configuration Error: /|\ Could not find Fisher matrix inverse after 10 attempts. Try with /o\ different fisher_tol, fisher_delta, input sigmas (or covmat), bestfit, or remove the option --fisher and run with Metropolis-Hastings or another sampling method.

I used this command: $ python montepython/MontePython.py run -o Planck18_50k_fisher -p input/base2018TTonly_lensing.param --bestfit bestfit/base2018TTonly_lensing.bestfit --method Fisher -N 50000 --superupdate 20 --conf mine.conf -f 1.5

without fisher, it worked fined.

[brinckmann]

Hi,

Thank you for reporting back.

The Fisher method isn't really very useful for data, as it computes the gradient away from the minimum and you don't know the true minimum of the likelihood, so I don't recommend it for that. It's more useful as a forecasting method to get an idea of what the sensitivity is like.

For dataset combinations where you don't have a starting covmat or good guesses for the sigmas to put in the param file you can sometimes try to start near the minimum with a very small jumping factor, e.g. -f 0.01 , and let it run for a while to try to accumulate enough points to force it to compute a not awful covmat with info and the --want-covmat flag and then start over in a new directory passing the new covmat with -c . Of course, that does require that you more or less know where the minimum is, as otherwise your chains will get stuck in a weird part of parameter space. Sometimes running a minimizer first can help with that, but they aren't always that great either unless you already more or less know the minimum. And it still might not work well if some of your guesses for the sigmas in your param file are way off. It's one of the potentially tricky problems of MCMC analysis.

Best, Thejs

[Kshitizpandey123]

Oh okay, I am actually going to use combinations now. This was very helpful.

Regards, Kshitiz

[brinckmann]

To clarify, if you have a good covmat for Planck alone and you add other data it's often enough to pass that covmat and run with the update and superupdate algorithms and let them adjust the covmat and jumping factor from your input guess.

But if that fails the suggestion I made with making the jumping factor small can sometimes help, if you know the best fitting area. Do be careful, because the jumping factor rescaling doesn't work well if you start with something small like 0.01, so after forcing it to compute a covmat after accumulating a lot of points you will usually want to start over in a new directory with a normal jumping factor (usually around 1.7-2.4, with 2.4 being for perfectly Gaussian posteriors).

Best, Thejs

[Kshitizpandey123]

Right, I got perfect Gaussian posteriors with a decent a.r. of 0.26. Thank you!

Now I am going to try to modify CLASS to analyze the Planck18 data with a dynamic Dark Energy Model.

[brinckmann]

Hi,

Great! Do make sure your maximum R-1 is about 0.01 or less as good rule of thumb for 4 or more chains (fewer and the R-1 value is artificially low) as the smoothing can trick a little and make you think things are converged before they are (but LCDM should be quick to get converged anyway). You can turn off smoothing with a plot file, there's an example in plot_files/example.plot listing three smoothing parameters (and more plotting options), which if you turn off you can see if the smoothed plots are a fair representation of the un-smoothed ones (which won't look as nice). This smoothing is common practice, but you don't want to misrepresent your data either. A big problem can be trying to smoothen non-Gaussian posteriors with a Gaussian fit, in which case you should turn off some of the smoothing parameters, but for LCDM they should be Gaussian.

You might need something different, but note that by default class has the CPL (w_0+w_a) parameterization, I think issues #32 , #212 and #280 is useful for how to set up. In short you need to set Omega_Lambda to 0 in order for the code to check for the fld parameters and it's better to sample w0_fld and the sum w0_fld+wa_fld, with appropriate bounds on the parameters, than the combination with wa_fld directly. If you need some other dynamical dark energy model that's, of course, something you'd rather need to look on the class side of things.

Best, Thejs