getplots gets stuck on MCMC output

[Carlos]

We managed to run a test of 300 points in our cluster using SBv1.35
with the MCMC option.
I tried to run getplots on the output, but it seems to get stuck:
cph@grad1> ./getplots GetPlotsSample.ini
Opening .info file: /home/cph/work/SBscans/output/run01.info
Number of parameters in MCMC chain: 71
Allocation
producing files with root run01/run01
Trying to read /home/cph/work/SBscans/output/run01.txt
Chain found, please wait...
outlier fraction 0.1433333
mean input multiplicity = 572.3900
single thin
using 300 rows, processing 4 parameters
All 2D plots status = F
effective number of samples (assuming indep): 31
Writing covariance matrix for 4 parameters
Best fit -Ln(like) = 7.67435
Mean(-Ln(like)) = mean(chisq/2) = 10.39384
-Ln(mean like) = 9.03288
Total number of variables = 6 Used variables = 4
Now doing 1D var j= 1 ix = 3 tails = 2 name = m_0 (GeV)
Opening for writing: plot_data/run01/run01_p1
Now doing 1D var j= 2 ix = 4 tails = 2 name = m_{1/2} (GeV)


And here it stays forever. Anything I am missing? The output files are
created.
Any hints?
Carlos.

[Roberto Trotta]

Hi Carlos

The issue is probably with the computation of the profile likelihood's
limits --- the code sometimes get stuck in an endless loop when trying
to figure out eg the 95% limits from the profile like if the
likelihood does not drop to below the 95% level within the points you
have. This is usually cured once you have more points (I will patch
this in a future release, it's very annoying, I agree).

In any case 300 points is way too small to get anything sensible. You
should aim for at list 3x10^4 (and discard the first couple of 1000 or
so for burnin) using MCMC (in our previous papers we actually used
3x10^5). We reccommend to use MultiNest instead, which is *much* more
efficient. Notice that if you do use MN you will not be able to run
getplots until the run has finished, because differently from MCMC
with MN the correctly weighted list of samples is only produced at the
end, once the sampler has evaluated the global evidence.

Roberto

[Carlos]

>
> In any case 300 points is way too small to get anything sensible.

yes I know, we are just setting up the programs and testing that we
can
compile and run the whole thing. No physics expected from such test.

> should aim for at list 3x10^4 (and discard the first couple of 1000 or
> so for burnin) using MCMC (in our previous papers we actually used
> 3x10^5). We reccommend to use MultiNest instead, which is *much* more
> efficient. Notice that if you do use MN you will not be able to run
> getplots until the run has finished, because differently from MCMC
> with MN the correctly weighted list of samples is only produced at the
> end, once the sampler has evaluated the global evidence.

ok, thanks for your answer.

Carlos.

[Farhan Feroz]

Hi Carlos,

Actually there is a workaround with MultiNest to produce the plots
using getplots. MultiNest produces the posterior output files after
every 1000 iterations which can be used with getplots. The results
would be a bit noisy in the initial stages of sampling but will get
better as the sampling proceeds. At the end of the sampling you would
get the actual posterior plots with getplots for the number of live
points used.

Farhan

[Roberto Trotta]

Hi

as a temporary quick-fix patch, I suggest the following:

in source/mcsamples.f90, replace the

function ConfidVal

with the code below. This will make sure the code does not get stuck
in an infinite loop (and will print out a warning in case the
convergence to the desired limit has not been attained).

I will fix this in a more sophisticated way in a future release with a
series of cumulative patches for smaller bugs like this that we have
found in the meantime.

Let me know if this helps,

Roberto



function ConfidVal(ix,limfrac,upper,ix1,ix2)
integer, intent(IN) :: ix
real, intent(IN) :: limfrac
logical, intent(IN) :: upper
integer, intent(IN), optional :: ix1,ix2
real ConfidVal, samps
real try, lasttry, try_t, try_b, try_t_prev, try_b_prev
integer l,t, jj
integer:: maxjj=1000
real*8, parameter :: tol = 1E-6, tolt=1.0
!find upper and lower bounds
jj =0
l=0
t=nrows-1
if (present(ix1)) l=ix1
if (present(ix2)) t = ix2
try_b = minval(coldata(ix,l:t))
try_t = maxval(coldata(ix,l:t))
samps = sum(coldata(1,l:t))

lasttry = -1
if (upper) then
do
try = sum(coldata(1,l:t),mask = coldata(ix,l:t) > (try_b +
try_t)/2)
if (try > samps*limfrac) then
try_b = (try_b+try_t)/2
else
try_t = (try_b+try_t)/2
end if
if (try == lasttry) exit
lasttry = try
end do
else
do
try = sum(coldata(1,l:t),mask = coldata(ix,l:t) < (try_b +
try_t)/2)
!write(*,'(A,E30.20)') 'Current mean value/target = ', try/
(samps*limfrac)
!write(*,'(A,E30.20)') 'This is for a mean = ', (try_b +
try_t)/2
if (try > samps*limfrac) then
try_t_prev = try_t
try_t = (try_b+try_t)/2
!write(*,'(A,E30.20)') 'Reducing try top = ', try_t
else
try_b_prev = try_b
try_b = (try_b+try_t)/2
if (try_b .eq. 0d0) then
try_b = try_b + 0.05*try_t
!write(*,*) 'Reset try_b to = ', try_b
end if
!write(*,'(A,E30.20)') 'Increasing try bot to = ', try_b
end if
!if ( (try == lasttry) .and. (abs(1 - try/(samps*limfrac))
< tol)) exit


!if ( (try == lasttry) .and. (abs(1 - try/(samps*limfrac))
> tol)) then
! write(*,*) 'sometjhing is rotten', (abs(1 - try/
(samps*limfrac)))
! write(*,*) 'latest: ', try_b, try_t
! write(*,*) 'previous: ', try_b_prev, try_t_prev
! stop
! exit
!end if
!if (abs(1 - try/(samps*limfrac)) < tol) then
! write(*,*) 'tolerance reacghed =', abs(1 - try/
(samps*limfrac))
! try_t = try
! exit
!end if

!write(*,*) jj, 'sometjhing is rotten? this has to go to
zero', (abs(1 - try/(samps*limfrac)))
!write(*,*) 'latest: ', try_b, try_t
!write(*,*) 'previous: ', try_b_prev, try_t_prev
!write(*,'(A,2E30.20)') '1E10*(try-last)%:', abs(1-try/
lasttry), tol
!write(*,'(A,2E30.20)') 'target :', abs(1 - try/
(samps*limfrac)), tolt
if (( abs(1-try/lasttry) < tol) .and. (abs(1 - try/
(samps*limfrac)) < tolt)) exit
lasttry = try
jj = jj+1
if(jj>maxjj) then
write(*,*) 'Warninng, could not converge to correct
limit! Continuing...'
exit
end if
end do
!write(*,*) 'DONE MINE LOWER JOB', try_t
!STOP
end if

ConfidVal = try_t

end function ConfidVal