Running CosmoMC to find the best fit point using only the supernovae data: Trail parameters excluded by prior or error in likelihood

[Minghua Li]

Hello,everyone!
I want to use the CosmoMC to find the best fit point (omegam, omegav) or (omegam, omegalambda). I didn't change anything in the params.ini in the cosmomc directory except the following two options:

1. change the "action" from 0 to 2, in order to find the best fit point only
2. make the program to use the supernovae data only, because I'm running it on a Duo Core E7200 PC. So I set :
use_CMB = F
use_HST = F
use_mpk = F
use_clusters = F
use_BBN = F
use_Age_Tophat_Prior = F
use_SN = T
use_lya = F
use_min_zre = 0

And no more else. When I compiled and ran the program I got the following output on the screen:

Random seeds: 5648, 18100 rand_inst: 0
Computing tensors: F
Doing CMB lensing: F
lmax = 2100
Number of C_ls = 3
num_params=13,index_nuisance=14
Varying 7 parameters ( 0 fast)
reading: params_CMB.covmat
Finding max-like point
params_used( 1)= 1
vector( 1)=
22.30000
params_used( 2)= 2
vector( 2)=
10.50000
params_used( 3)= 3
vector( 3)=
519.9999
params_used( 4)= 4
vector( 4)=
3.000000
params_used( 5)= 8
vector( 5)=
95.00000
params_used( 6)=11
vector( 6)=
300.0000
params_used( 7)=13
vector( 7)=
2.500000
Reading: supernovae data
ffn=
155.8963
ffn=
155.8963
ffn=
155.9300
ffn=
155.8628
ffn=
155.9639
ffn=
155.8294
ffn=
156.0333
...

ffn=
155.8963
ffn=
155.8963
ffn=
155.8963
ffn=
155.8963
ffn=
155.8963
ffn=
155.8963
Print the value of Scales%PMax(i) , Params%P(i) and Scales%PMin(i) in sequence , with i from 1 to 7:
1 0.1000000
1 2.2299999E-02
1 4.9999999E-03
2 0.9900000
2 0.1050000
2 9.9999998E-03
3 10.00000
3 1.040000
3 0.5000000
4 0.8000000
4 -5.9999940E-03
4 9.9999998E-03
5 1.500000
5 0.9500000
5 0.5000000
6 4.000000
6 3.000000
6 2.700000
7 2.000000
7 1.000000
7 0.0000000E+00

Attention: there are any(Params%P > Scales%PMax) .or. any(Params%P < Scales%PMin!!
About to set GetLogLike = logZero!!!
ffn=
1.0000000E+30
ERROR: Trial parameters excluded by prior or error in likelihood
Try starting further away from problem regions?

I have some questions:
1. What is the sn_covmat_sys.txt in the supernovae.f90 ?
Is it the covariance matrix or something? If it is ,could you please explain more about how it is related to the sn_moduli ? Because I couldn't find any useful information anywhere in any of my statistics text book or even Wiki encyclopedia online.
2. What is the real array P , which is defined in the Type ParamSet in CMB_Cls_simple.f90 , used for?
Is it used to give each parameter a number so it would be much more easy to manipulate it in the conjgrad_wrapper.f90 using the "multiple dimension conjurgate gradient method" ? Because it seems to make the parameters much difficult to trace!!
3. Last but not least, where does the abnormal Params%P(4) (See the error message above) stem from?
Because I tried using another set of supernovae data, and it gave different values of ffn, but the same abnormal values of Scales%PMax(i) , Params%P(i) and Scales%PMin(i) . That makes me believe that the error doesn't lie in the supernovae data----it stems from the program, it stems from trying to assign a number to the real P array in the Type ParamSet !!! That brings up the Question 2.


Could anyone please help me out of this, if you're patient enough to finished reading this tedious error report? I will continue trying myself while waiting. I would add newpost when I know more,
Thanks for your attention !!! Best regards, limh

[Minghua Li]

I found something that may be useful :

Inline =Ini_Read_String (numcat('param', i ), .true. )
Read(Inline, *, err= 100) center, Scales%PMin(i), Scales%PMax(i), wid, Scales%PWidth(i)

The above two lines are in the paramdef.f90 . It was meant to read in the values of center, Scales%PMin(i), etc. which, I believe, are related to the starting point of the Maximum Likelihood Estimation (MLE) Procedure. [ Model fitting is an iterative procedure: the user has to specify a set of starting values for the parameters (essentially an initial 'first guess') which the optimisation algorithm will take and try to improve on. ]
Moreover, see below( also adapted from paramdef.f90 ):

if (new_chains) then ! new_chains was set to be .true.
do ! Try to get acceptable values in range
if (wid < 0) then
!This case we want half gaussian, width -wid
!e.g. for positive definite parameters
Params%P(i) = center - abs(Gaussian1())*wid
else
Params%P(i) = center + Gaussian1()*wid
end if
!Repeat until get acceptable values in range
if (Params%P(i)>= Scales%PMin(i) .and. Params%P(i) <= Scales%PMax(i)) exit
end do

end if
end do


The point is : the above codes has already get the Params%P(i) in the range [Scales%PMin(i) , Scales%PMax(i) ], why would it later aquire an abnormal value anb be out of this range? [See the output in the first post.]

The problem lies in the funtion ffn in conjgrad_wrapper.f90, it did some modifications to P%P(ii) , which is just Params%P(i):

function ffn(vect,OtherParams, keepNew)
use ParamDef
use CalcLike
implicit none
logical, intent(in) :: KeepNew

Type(ParamSet) P, OtherParams
integer i,ii
real ffn
real vect(num_params_used)

! OtherParams contains all the other info about the cosmology
! The subroutines only know about vect, and they make changes to this.
! But they carry around OtherParams (never changing it) so it can be used here.

! scale the params so they are all roughly the same order of magnitude
P = OtherParams

do i=1, num_params_used
ii=params_used(i)
P%P(ii)=vect(i)*Scales%PWidth(ii)
end do
...

It was this part of codes that generated an abnormal Params%P(i) !

But why would it be?! Is there anything wrong with the procedure, because I didn't change anything except for turning on the "Use_SN" (Use_SN= T) and turning off all other "Use_xxx"(= F), and setting "action" to be 2 ?

All help is welcome, thanks!!

[Minghua Li]

I found something new, see the feedback of my program:

....

Call dffn..........................................................

Finding the gradient numerically...
Print out the vector(1:7), Scales%PWidth(1:7 )and P%P(1:7) for check:
vector( 1)=
22.30000
vector( 2)=
10.50000
vector( 3)=
519.9999
vector( 4)=
3.000000
vector( 5)=
95.00000
vector( 6)=
300.0000
vector( 7)=
2.500000
Scales%PWidth(1:7)=
1 1.0000000E-03
2 9.9999998E-03
3 2.0000001E-03
4 2.9999999E-02
5 9.9999998E-03
6 9.9999998E-03
7 0.4000000
P%P(1:7)=
1 2.2299999E-02
2 0.1050000
3 1.040000
4 9.0000004E-02
5 0.9500000
6 3.000000
7 1.000000
Way 2!
SN chisq: 311.7926
funcvect=ffn=
155.8963
*******************************************************************************
**********
The 1 round of do loop
Print out the vector(1:7), Scales%PWidth(1:7 )and P%P(1:7) for check:
vector( 1)=
22.50000
vector( 2)=
10.50000
vector( 3)=
519.9999
vector( 4)=
3.000000
vector( 5)=
95.00000
vector( 6)=
300.0000
vector( 7)=
2.500000
Scales%PWidth(1:7)=
1 1.0000000E-03
2 9.9999998E-03
3 2.0000001E-03
4 2.9999999E-02
5 9.9999998E-03
6 9.9999998E-03
7 0.4000000
P%P(1:7)=
1 2.2500001E-02
2 0.1050000
3 1.040000
4 9.0000004E-02
5 0.9500000
6 3.000000
7 1.000000
Way 2!
SN chisq: 311.8600
Print out the vector(1:7), Scales%PWidth(1:7 )and P%P(1:7) for check:
vector( 1)=
22.10000
vector( 2)=
10.50000
vector( 3)=
519.9999
vector( 4)=
3.000000
vector( 5)=
95.00000
vector( 6)=
300.0000
vector( 7)=
2.500000
Scales%PWidth(1:7)=
1 1.0000000E-03
2 9.9999998E-03
3 2.0000001E-03
4 2.9999999E-02
5 9.9999998E-03
6 9.9999998E-03
7 0.4000000
P%P(1:7)=
1 2.2100000E-02
2 0.1050000
3 1.040000
4 9.0000004E-02
5 0.9500000
6 3.000000
7 1.000000
Way 2!
SN chisq: 311.7255

....
SN chisq: 310.0932
gradient wrt 1 th element is: 0.1641214
*******************************************************************************
**********
....


The 4 round of do loop
FORTRAN PAUSE
PAUSE prompt>
Print out the vector(1:7), Scales%PWidth(1:7 )and P%P(1:7) for check:
vector( 1)=
22.30000
vector( 2)=
10.50000
vector( 3)=
519.9999
vector( 4)=
3.200000
vector( 5)=
95.00000
vector( 6)=
300.0000
vector( 7)=
2.500000
Scales%PWidth(1:7)=
1 1.0000000E-03
2 9.9999998E-03
3 2.0000001E-03
4 2.9999999E-02
5 9.9999998E-03
6 9.9999998E-03
7 0.4000000
P%P(1:7)=
1 2.2299999E-02
2 0.1050000
3 1.040000
4 9.6000008E-02
5 0.9500000
6 3.000000
7 1.000000
Way 2!
SN chisq: 311.7926
...


Print out the vector(1:7), Scales%PWidth(1:7 )and P%P(1:7) for check:
vector( 1)=
22.30000
vector( 2)=
10.50000
vector( 3)=
519.9999
vector( 4)=
-0.1999998
vector( 5)=
95.00000
vector( 6)=
300.0000
vector( 7)=
2.500000
Scales%PWidth(1:7)=
1 1.0000000E-03
2 9.9999998E-03
3 2.0000001E-03
4 2.9999999E-02
5 9.9999998E-03
6 9.9999998E-03
7 0.4000000
P%P(1:7)=
1 2.2299999E-02
2 0.1050000
3 1.040000
4 -5.9999940E-03
5 0.9500000
6 3.000000
7 1.000000

1 0.1000000
1 2.2299999E-02
1 4.9999999E-03
2 0.9900000
2 0.1050000
2 9.9999998E-03
3 10.00000
3 1.040000
3 0.5000000
4 0.8000000
4 -5.9999940E-03
4 9.9999998E-03
5 1.500000
5 0.9500000
5 0.5000000
6 4.000000
6 3.000000
6 2.700000
7 2.000000
7 1.000000
7 0.0000000E+00
There are any(Params%P > Scales%PMax) .or. any(Params%P < Scales%PMin, about to
set GetLogLike = logZero!!!
ffn=
1.0000000E+30
About to call outofbound in ffn!!!
FORTRAN PAUSE
PAUSE prompt>
ERROR: Trial parameters excluded by prior or error in likelihood
Try starting further away from problem regions?



Notice that :
The value of vect(i) are swaying left and right per delta_vect(=0.2) in the i th round of do loop. It was the shifting of vect(i) that caused the shifting Params%P(i), which turned out to be out of [ Scales%PMin(i), Scales%PMax(i) ]!!!!

The question is : why? Is this error due to the improper initial value of the starting point in the parameters space? Or do I need more constraints and condition to rule out the "abnormal" vect(i) ?

[Antony Lewis]

As of the Oct 2012 version, best-fit finding (action=2) and estimate_propose_matrix should work more reliably, esp. with parameters that have hard prior cuts.