florianhartig / BayesianTools

General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics
https://cran.r-project.org/web/packages/BayesianTools/index.html
117 stars 29 forks source link

Start values be outside prior range in DREAMzs / DEzs #189

Open evodevosys opened 4 years ago

evodevosys commented 4 years ago

I'm repeatedly getting these two warnings:

4: Start values outside prior range
5: Z matrix values outside prior range

What does this mean? It seems to happen repeatedly during the run, so I don't think it's just when I start a chain. Is it something to worry about? Should I modify my priors or runMCMC call to fix it?

I've pasted my priors and call below:

Generic prior density and prior sampler calls. The currentParInfo global variable represents a list with information that is sent into the actual distribution functions. The funcs field is one of the distribution functions below.

prdensities=function(pars){
  ds=sum(sapply(1:length(pars),function(i){
    currentParInfo$funcs[[i]](a=currentParInfo$low[i],b=currentParInfo$high[i],m=currentParInfo$best[i],s=currentParInfo$s[i],p=pars[i])
  }))
  return(ds)
}

prsamplers=function(n=1){
  ds=sapply(1:length(currentParInfo$low),function(i){
    currentParInfo$funcs[[i]](a=currentParInfo$low[i],b=currentParInfo$high[i],m=currentParInfo$best[i],s=currentParInfo$s[i],n=n)
  })
  return(ds)
}

The distribution functions are below. Variables a and b are lower and upper bounds where the distribution requires truncation. These are the same values that are passed into createPrior(). m is the mean and is the standard deviation where the distribution requires it. Again, the same as passed into createPrior() [best=mean]

The beta distribution goes between 0 and 1. Best is passed in as 1. Could this be the problem?

# dtrunc and rtrunc are from the LaplacesDemon package
dist.normtrunc<-function(a,b,m,s,p=NULL,n=NULL){
  if(is.null(n)){
    return(dtrunc(p,'norm',a=a,b=b,mean=m,sd=s,log=TRUE))
  } else {
    return(rtrunc(n,'norm',a=a,b=b,mean=m,sd=s))
  }
}

dist.unif<-function(a,b,...,p=NULL,n=NULL){# the ... lets it take m and s and not throw and error even though it doesn't use them
  if(is.null(n)){
    return(dunif(p,min=a,max=b,log=TRUE))
  } else {
    return(runif(n,min=a,max=b))
  }
}

dist.beta<-function(...,p=NULL,n=NULL){# the ... lets it take a,b,m and s and not throw and error even though it doesn't use them
  # This is fixed at the beta(4,1) distribution
  if(is.null(n)){
    return(dbeta(p,shape1=4,shape2=1,log=TRUE))
  } else {
    return(rbeta(n,shape1=3,shape2=1))
  }
}

The runMCMC call is:

priors=createPrior(density=prdensities,sampler=prsamplers,lower=currentParInfo$low,upper=currentParInfo$high,best=currentParInfo$best)
  bSetup=createBayesianSetup(likelihood=runTheModelAndGetLogLikelihoods,prior=priors)
  settings.1=list(iterations=50000,burnin=10000,nrChains=3,adaptation=9000)#for DREAM

  results=runMCMC(bayesianSetup =bSetup,settings=settings.1,sampler='DREAMzs')
florianhartig commented 4 years ago

The most likely reason is that your prior sampler produces values that are not within your prior or within the lower / upper conditions that you specify in the BayesianSetup.

In particular the lower upper has been an issue in the past - if this is inconsistent with your prior, the z-Matrix could be filled with values that have no prior support.

evodevosys commented 4 years ago

Thanks for the quick reply!! While you were replying I edited my question that gives a little more detail about just that. The upper and lower bounds for the sampler are truncated at the same values of low and high that I pass into createPrior(). [note that I don't use dist.norm() I'll take that out. The uniform, beta, and truncated normal are all bounded]

florianhartig commented 4 years ago

You're sure? I'd bet your sampler creates samples outside the prior range. Try x = bSetup$prior$sampler(1000), and then bSetup$prior$density(x) and check if there are any -Inf values.

evodevosys commented 4 years ago

I tried it with 100000 values. No -Inf

> x=bSetup$prior$sampler(100000)
> any(is.infinite(bSetup$prior$density(x)))
[1] FALSE
florianhartig commented 4 years ago

That is weird. The message is quite clear, and it all sounds as if your prior sampler produces values that are outside the prior range.

If you look into the Z-matrix of your sampler (directly after starting it up), e.g. like that

settings = list(iterations = 4)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DREAMzs", settings = settings)
out$Z

what kind of values do you see? Are they really inside the prior range? Do a visual check of prior and likelihood of these values - anything unusual?

If nothing of this helps, I guess the only thing is to look into the code and see what's going wrong. If you can send me a reproducible example of what you're doing, I could have a look, otherwise switch the debug mode on check why this error is occurring.

To be clear - it's warning, in principle the algorithm should still work, but it's a weird thing that this happens, so I would investigate what's going on.

evodevosys commented 4 years ago

The out$Z values are all inside the prior range. The likelihoods and priors don't look particularly odd to me - but I haven't looked at a ton of these: image

image

I'll do some more investigating, see if I can narrow down exactly when it is happening and if there is anything consistent about the parameter values it is trying at those times, and put together a reproducible example

Thanks

florianhartig commented 4 years ago

Hi, yes, this looks completely normal.

I have no idea why you get this error, and I can't reproduce it with an example model. To investigate further, I would need a reproducible example.