CosmoCoffee Forum Index CosmoCoffee

 
 FAQFAQ   SearchSearch  MemberlistSmartFeed   MemberlistMemberlist    RegisterRegister 
   ProfileProfile   Log inLog in 
Arxiv New Filter | Bookmarks & clubs | Arxiv ref/author:

CosmoMC: problem with mock data
 
Post new topic   Reply to topic    CosmoCoffee Forum Index -> Computers and software
View previous topic :: View next topic  
Author Message
Paolo Campeti



Joined: 09 Jan 2018
Posts: 5
Affiliation: SISSA, Trieste

PostPosted: February 01 2018  Reply with quote

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:

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:
dataCols = 'TT EE TE'
fwhm_arcmin = 7.0
NoiseVar = 3e-4

and create the file MyForecast.ini
Code:
cmb_dataset[MyForecast]=data/MyForecast/simulatedCls_exactsim.dataset

I run CosmoMC with the test.ini file
Code:
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:
use_nonlinear_lensing = F
block_semi_fast = T
CMB_lensing = F
stop_on_error=  F

while the params_CMB_defaults.ini reads
Code:
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:
Number of chains used =  4
 var(mean)/mean(var), remaining chains, worst e-value: R-1 =       0.00241
RL: Thin for Markov:  27
RL: Thin for indep samples:   31
RL: Estimated burn in steps:  130  ( 52  rows)
using 55299 rows, 47 parameters; mean weight 2.52382502396, tot weight 139565.0
Approx indep samples (N/corr length): 4362.0
Equiv number of single samples (sum w)/max(w): 5169.0
Effective number of weighted samples (sum w)^2/sum(w^2): 32924
Best fit sample -log(Like) = 268.629400
Ln(mean 1/like) = 276.122745
mean(-Ln(like)) = 271.872252
-Ln(mean like)  = 270.858158

and the chains have around 30000 lines. Apparently the chains remain stuck, as can be seen from the triangle plots below

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.
Back to top
View user's profile  
Antony Lewis



Joined: 23 Sep 2004
Posts: 1342
Affiliation: University of Sussex

PostPosted: February 02 2018  Reply with quote

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.
Back to top
View user's profile [ Hidden ] Visit poster's website
Paolo Campeti



Joined: 09 Jan 2018
Posts: 5
Affiliation: SISSA, Trieste

PostPosted: February 02 2018  Reply with quote

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:
 Number of MPI processes:           4
 file_root:test_the_likelihood
 Random seeds:  1630,  8395 rand_inst:   1
 Random seeds:  1937,  8396 rand_inst:   4
 Random seeds:  1733,  8395 rand_inst:   2
 Random seeds:  1841,  8396 rand_inst:   3
 Doing non-linear Pk: F
 Doing CMB lensing: F
 TT lmax =  2500
 EE lmax =  2500
 ET lmax =  2500
 lmax_computed_cl  =  2500
 Computing tensors: F
 max_eta_k         =    6250.000   
 transfer kmax     =    1.000000   
 adding parameters for: simulatedCls_exactsim
 Fast divided into            1  blocks
  6 parameters ( 6 slow ( 2 semi-slow),  0 fast ( 0 semi-fast))
 skipped unused params: calPlanck acib217 xi asz143 aps100 aps143 aps143217 aps217 aksz kgal100 kgal143 kgal143217 kgal217 cal0 cal2
 skipped unused params: calPlanck acib217 xi asz143 aps100 aps143 aps143217 aps217 aksz kgal100 kgal143 kgal143217 kgal217 cal0 cal2
 skipped unused params: calPlanck acib217 xi asz143 aps100 aps143 aps143217 aps217 aksz kgal100 kgal143 kgal143217 kgal217 cal0 cal2
 skipped unused params: calPlanck acib217 xi asz143 aps100 aps143 aps143217 aps217 aksz kgal100 kgal143 kgal143217 kgal217 cal0 cal2
        Time for theory:    2.87372
 Time for simulatedCls_exactsim:   1.126313209533691E-002
    loglike     chi-sq
  20121.235  40242.470   CMB: MyForecast = simulatedCls_exactsim
        Time for theory:    2.89992
 Time for simulatedCls_exactsim:   9.938955307006836E-003
    loglike     chi-sq
  20121.235  40242.470   CMB: MyForecast = simulatedCls_exactsim
Test likelihoods done, total logLike, chi-eq =   20121.235  40242.470
 Likelihood calculation time (seconds)= 2.9046
Test likelihoods done, total logLike, chi-eq =   20121.235  40242.470
 Likelihood calculation time (seconds)= 2.9294
        Time for theory:    2.96636
        Time for theory:    2.96692
 Time for simulatedCls_exactsim:   1.102399826049805E-002
    loglike     chi-sq
  20121.235  40242.470   CMB: MyForecast = simulatedCls_exactsim
 Time for simulatedCls_exactsim:   9.715080261230469E-003
    loglike     chi-sq
  20121.235  40242.470   CMB: MyForecast = simulatedCls_exactsim
Test likelihoods done, total logLike, chi-eq =   20121.235  40242.470
 Likelihood calculation time (seconds)= 2.9970
Test likelihoods done, total logLike, chi-eq =   20121.235  40242.470
 Likelihood calculation time (seconds)= 2.9963
Total time:  3  (   0.00086 hours  )

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
Back to top
View user's profile  
Antony Lewis



Joined: 23 Sep 2004
Posts: 1342
Affiliation: University of Sussex

PostPosted: February 05 2018  Reply with quote

Looks like your parameters are quite difference. Check the .inputparams file actually lists all as you expected for the true simulated model.
Back to top
View user's profile [ Hidden ] Visit poster's website
Paolo Campeti



Joined: 09 Jan 2018
Posts: 5
Affiliation: SISSA, Trieste

PostPosted: February 05 2018  Reply with quote

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:
dataCols = 'TT TE EE BB PP'

to
Code:
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:
fwhm_arcmin = 5.
NoiseVar = 4e-5
ENoiseFac = 4
lmin = 2
lmax = 2500
fsky = 0.57

to
Code:
fwhm_arcmin = 7.0
NoiseVar = 3e-4
ENoiseFac = 2
lmin = 2
lmax = 2000
fsky = 1.0
)
Back to top
View user's profile  
Antony Lewis



Joined: 23 Sep 2004
Posts: 1342
Affiliation: University of Sussex

PostPosted: February 06 2018  Reply with quote

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.
Back to top
View user's profile [ Hidden ] Visit poster's website
Display posts from previous:   
Post new topic   Reply to topic    CosmoCoffee Forum Index -> Computers and software All times are GMT + 5 Hours
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © 2001, 2005 phpBB Group. Sponsored by WordWeb online dictionary and dictionary software.