CosmoMC: problem with mock data

Use of Cobaya. camb, CLASS, cosmomc, compilers, etc.
Post Reply
Paolo Campeti
Posts: 5
Joined: January 09 2018
Affiliation: SISSA, Trieste

CosmoMC: problem with mock data

Post by Paolo Campeti » February 01 2018

Hello everyone,
I'm trying to use CosmoMC-Nov2016 with only simulated Planck data generated using the makePerfectForecastDataset.py script but the parameters recovered through MCMC are very different from the ones used to generate the data. I use CAMB (the version included into CosmoMC-Nov2016) to produce the scalar Cls, modifying only the following parameters in the params.ini file (all the other parameters are the default ones with which CosmoMC is dowloaded)

Code: Select all

l_max_scalar = 2000
do_lensing = F
use_physical = T
ombh2 = 0.022
omch2 = 0.11
omnuh2 = 0.0
omk = 0
hubble = 70
re_optical_depth = 0.09
scalar_spectral_index(1) = 0.96
scalar_amp(1) = 2.4e-9
l_sample_boost = 50
Then I modify makePerfectForecastDataset.py in the following lines

Code: Select all

dataCols = 'TT EE TE'
fwhm_arcmin = 7.0
NoiseVar = 3e-4
and create the file MyForecast.ini

Code: Select all

cmb_dataset[MyForecast]=data/MyForecast/simulatedCls_exactsim.dataset
I run CosmoMC with the test.ini file

Code: Select all

DEFAULT(batch2/MyForecast.ini)
DEFAULT(batch2/common.ini)
MPI_Max_R_ProposeUpdate = 30
propose_matrix= planck_covmats/base_TT_lowTEB_plik.covmat
action = 0
start_at_bestfit =F
feedback=1
use_fast_slow = T
checkpoint = T
sampling_method = 7
dragging_steps  = 3
propose_scale = 2
indep_sample=10
My common.ini file is only modified with

Code: Select all

use_nonlinear_lensing = F
block_semi_fast = T
CMB_lensing = F
stop_on_error=  F
while the params_CMB_defaults.ini reads

Code: Select all

param[omegabh2] = 0.0221 0.005 0.1 0.0001 0.0001
param[omegach2] = 0.12 0.01 0.2 0.001 0.0005
param[theta] = 1.03 0.5 1.2 0.0004 0.0002
param[tau] = 0.089 0.01 0.3 0.01 0.005
neutrino_hierarchy = degenerate
num_massive_neutrinos=1
param[mnu] = 0.0
param[meffsterile] = 0
param[omegak] = 0
param[w] = -1
param[nrun] = 0 
param[nrunrun] = 0
param[r] = 0
param[wa] = 0
param[nnu] = 0.0
param[yhe] = 0.24
param[alpha1] = 0
param[deltazrei] = 0.5
param[Alens] = 1
param[Alensf]=-1
param[fdm] = 0
param[ns] = 0.96 0.8 1.2 0.004 0.002
param[logA] = 3.1 2 4 0.001 0.001
param[Aphiphi] = 1 
After convergence (R-1<0.01) In the .margestats file the mean values are very far from the fiducial ones used to generate the mock data, for instance
omegabh2 ~ 0.019 insted of 0.022
omegach2 ~ 0.069 instead of 0.11
theta ~ 1.04
tau~0.0757 instead of 0.089
logA~2.94 instead of 3.13
ns~0.801 instead of 0.96
H0*~46.7 instead of 70
omegal*~0.592 instead of 0.73
GetDist gives

Code: Select all

Number of chains used =  4
 var&#40;mean&#41;/mean&#40;var&#41;, remaining chains, worst e-value&#58; R-1 =       0.00241
RL&#58; Thin for Markov&#58;  27
RL&#58; Thin for indep samples&#58;   31
RL&#58; Estimated burn in steps&#58;  130  &#40; 52  rows&#41;
using 55299 rows, 47 parameters; mean weight 2.52382502396, tot weight 139565.0
Approx indep samples &#40;N/corr length&#41;&#58; 4362.0
Equiv number of single samples &#40;sum w&#41;/max&#40;w&#41;&#58; 5169.0
Effective number of weighted samples &#40;sum w&#41;^2/sum&#40;w^2&#41;&#58; 32924
Best fit sample -log&#40;Like&#41; = 268.629400
Ln&#40;mean 1/like&#41; = 276.122745
mean&#40;-Ln&#40;like&#41;&#41; = 271.872252
-Ln&#40;mean like&#41;  = 270.858158
and the chains have around 30000 lines. Apparently the chains remain stuck, as can be seen from the triangle plots below
Image
I tried also increasing the temperature to 20 but I obtained the same wrong parameters. Should I add a prior on H0 maybe? Any help is greatly appreciated, many thanks for your time and patience.

Antony Lewis
Posts: 1941
Joined: September 23 2004
Affiliation: University of Sussex
Contact:

Re: CosmoMC: problem with mock data

Post by Antony Lewis » February 02 2018

Are you plotting the burn in period there? Make sure getdist is set to ingore about the first 30% of the chains.

Not sure what the bias is, though I've not tested unlensed for years. Use action=4, with parameters centred at the true values, to see if you get out a chi2 of zero, and set test_output_root to get a dump of what it is calculating.

Paolo Campeti
Posts: 5
Joined: January 09 2018
Affiliation: SISSA, Trieste

CosmoMC: problem with mock data

Post by Paolo Campeti » February 02 2018

Thank you very much for your answer Prof.Lewis,
No in the plots above ignore_rows was already set to 0.3 in the distparams.ini file.
Using action=4 with parameters set exactly to the true values I obtain the following output:

Code: Select all

 Number of MPI processes&#58;           4
 file_root&#58;test_the_likelihood
 Random seeds&#58;  1630,  8395 rand_inst&#58;   1
 Random seeds&#58;  1937,  8396 rand_inst&#58;   4
 Random seeds&#58;  1733,  8395 rand_inst&#58;   2
 Random seeds&#58;  1841,  8396 rand_inst&#58;   3
 Doing non-linear Pk&#58; F
 Doing CMB lensing&#58; F
 TT lmax =  2500
 EE lmax =  2500
 ET lmax =  2500
 lmax_computed_cl  =  2500
 Computing tensors&#58; F
 max_eta_k         =    6250.000    
 transfer kmax     =    1.000000    
 adding parameters for&#58; simulatedCls_exactsim
 Fast divided into            1  blocks
  6 parameters &#40; 6 slow &#40; 2 semi-slow&#41;,  0 fast &#40; 0 semi-fast&#41;&#41;
 skipped unused params&#58; calPlanck acib217 xi asz143 aps100 aps143 aps143217 aps217 aksz kgal100 kgal143 kgal143217 kgal217 cal0 cal2
 skipped unused params&#58; calPlanck acib217 xi asz143 aps100 aps143 aps143217 aps217 aksz kgal100 kgal143 kgal143217 kgal217 cal0 cal2
 skipped unused params&#58; calPlanck acib217 xi asz143 aps100 aps143 aps143217 aps217 aksz kgal100 kgal143 kgal143217 kgal217 cal0 cal2
 skipped unused params&#58; calPlanck acib217 xi asz143 aps100 aps143 aps143217 aps217 aksz kgal100 kgal143 kgal143217 kgal217 cal0 cal2
        Time for theory&#58;    2.87372
 Time for simulatedCls_exactsim&#58;   1.126313209533691E-002
    loglike     chi-sq
  20121.235  40242.470   CMB&#58; MyForecast = simulatedCls_exactsim
        Time for theory&#58;    2.89992
 Time for simulatedCls_exactsim&#58;   9.938955307006836E-003
    loglike     chi-sq
  20121.235  40242.470   CMB&#58; MyForecast = simulatedCls_exactsim
Test likelihoods done, total logLike, chi-eq =   20121.235  40242.470
 Likelihood calculation time &#40;seconds&#41;= 2.9046
Test likelihoods done, total logLike, chi-eq =   20121.235  40242.470
 Likelihood calculation time &#40;seconds&#41;= 2.9294
        Time for theory&#58;    2.96636
        Time for theory&#58;    2.96692
 Time for simulatedCls_exactsim&#58;   1.102399826049805E-002
    loglike     chi-sq
  20121.235  40242.470   CMB&#58; MyForecast = simulatedCls_exactsim
 Time for simulatedCls_exactsim&#58;   9.715080261230469E-003
    loglike     chi-sq
  20121.235  40242.470   CMB&#58; MyForecast = simulatedCls_exactsim
Test likelihoods done, total logLike, chi-eq =   20121.235  40242.470
 Likelihood calculation time &#40;seconds&#41;= 2.9970
Test likelihoods done, total logLike, chi-eq =   20121.235  40242.470
 Likelihood calculation time &#40;seconds&#41;= 2.9963
Total time&#58;  3  &#40;   0.00086 hours  &#41;
a very high chi square.
I tried also to plot the temperature simulated Cls (in blue) together with the test_output_root .theory_cl Cls (in yellow):
Clsexactsim_and_output_theory.cl.png
also the TE Cls
TEsimandtheoryCls.png
and the EE Cls
EEsimandtheory.png

Antony Lewis
Posts: 1941
Joined: September 23 2004
Affiliation: University of Sussex
Contact:

Re: CosmoMC: problem with mock data

Post by Antony Lewis » February 05 2018

Looks like your parameters are quite difference. Check the .inputparams file actually lists all as you expected for the true simulated model.

Paolo Campeti
Posts: 5
Joined: January 09 2018
Affiliation: SISSA, Trieste

CosmoMC: problem with mock data

Post by Paolo Campeti » February 05 2018

Thank you very much Prof. Lewis for your time and patience,
The .inputparams file
test_the_likelihood.inputparams
looks ok as all the parameters here have the values I fixed in the .ini files.
Is it possible that I'm making some error in the way I modified the makePerfectForecastDataset.py script?
My input Cls are in the form of 3 columns: TT EE TE so I changed the line

Code: Select all

dataCols = 'TT TE EE BB PP'
to

Code: Select all

dataCols = 'TT EE TE'
Is this correct? Maybe the columns order gives some problem where CAMB is called by CosmoMC to compute the test likelihood?
(I also modified the following values in makePerfectForecastDataset.py from

Code: Select all

fwhm_arcmin = 5.
NoiseVar = 4e-5
ENoiseFac = 4
lmin = 2
lmax = 2500
fsky = 0.57

to

Code: Select all

fwhm_arcmin = 7.0
NoiseVar = 3e-4
ENoiseFac = 2
lmin = 2
lmax = 2000
fsky = 1.0
)

Antony Lewis
Posts: 1941
Joined: September 23 2004
Affiliation: University of Sussex
Contact:

Re: CosmoMC: problem with mock data

Post by Antony Lewis » February 06 2018

it doesn't look to me as though your .inputparam central (first) values for the param[xx] settings agree with what you simulated, which is what you need for action=4 to check out. Esp. check the theta values are consistent as well.

Post Reply