Truncated normal distribution with WinBUGS

Question

I am trying to do a very simple regression analysis using a truncated normal distribution with WinBUGS.

Let's suppose the following model $Y|\beta_0,\beta_1,\sigma \sim Normal(\beta_0+\beta_1X,\sigma)$ where $Y$ is the grade of 13 students (from 0 to 10) and $X$ is the gender of each student (1 for Males and 0 for Females). We wonder whether the grade depends on the gender, so we could state if males have higher grades than females or the other way around. The following code is fine and does the job. However, as the grade is scaled from 0 to 10, I thought that using a truncated normal distribution could improve the model: $Y|\beta_0,\beta_1,\sigma \sim TrNormal(\beta_0+\beta_1X,\sigma,\alpha=0,\beta=10)$ where $\alpha$ and $\beta$ are the lower and upper bounds respectively.

The truncated normal distribution for WinBUGS is available as a WBDev shared component. However, when I am trying to use it with the following script it hangs unexpectedly. It returns a messages titled "undefined real result". What am I doing wrong?

WinBUGS script:

#MODEL
model
{
    for(i in 1:N)
    {
       mu[i]<-b0+b1*x[i]
       #y[i]~dnorm(mu[i],tau)
       y[i]~djl.dnorm.trunc(mu[i],tau,0,10)
    }
    tau~dgamma(0.01,0.01)
    b0~dnorm(0,1.0E-6)
    b1~dnorm(0,1.0E-6)
    sigma<-sqrt(1/tau)
}
#DATA
list(N=13,y=c(9.6,7.0,5.0,8.0,8.4,6.4,6.1,9.1,8.8,5.7,8.9,6.1,6.5),x=c(1,1,1,1,1,1,0,0,0,0,0,0,0))
#INITS
list(tau=1, b0=0, b1=0)

Answer 1

In openbugs (so I would assume the same in winbugs, but could be wrong) you can do a truncated distribution using the T function, something like:

y[i] ~ dnorm(x[i],tau)T(0,10)

Answer 2

I found the solution by consulting the WinBUGS User Manual (Page 47). I had to modify the initial parameters but I end up using the I() function with dnorm.

undefined real result indicates numerical overflow. Possible reasons include:

• initial values generated from a 'vague' prior distribution may be numerically extreme - specify appropriate initial values;

• numerically impossible values such as log of a non-positive number - check, for example, that no zero expectations have been given when Poisson modelling;

• numerical difficulties in sampling. Possible solutions include:

  • better initial values;
  • more informative priors - uniform priors might still be used but with their range restricted to plausible values;
  • better parameterisation to improve orthogonality;
  • standardisation of covariates to have mean 0 and standard deviation 1.
• can happen if all initial values are equal.

Probit models are particularly susceptible to this problem, i.e. generating undefined real results. If a probit is a stochastic node, it may help to put reasonable bounds on its distribution, e.g.

probit(p[i]) <- delta[i] delta[i] ~ dnorm(mu[i], tau)I(-5, 5)

This trap can sometimes be escaped from by simply clicking on the update button. The equivalent construction

p[i] <- phi(delta[i])

may be more forgiving.

Answer 3

According to http://www.unc.edu/courses/2010fall/ecol/563/001/docs/solutions/assign10.htm, two gotchas that may apply are that

1. the upper truncation point may need to be adjusted (e.g., from 10 to 10.0001)
2. initial values that seed the chains need to be specified explicitly, because WinBUGS may pick ones that lie outside of ($\alpha$,$\beta$).