Normal vs Truncated Normal

[fbud...@gmail.com]

Hi all,

We have been working on a model where our response variable can't be negative. We were using a normal distribution in the model development stage, but a truncated normal is a better match for the support our data. 

When we fit the non-truncated normal, we get an estimate of the mean that makes sense given the data. However, when we fit the truncated normal, the posterior samples for the mean are all pushed up against zero, the mean doesn't represent the mean of the data, and in more complex cases this produces some very abnormal posterior distributions.

Maybe this is a feature of the truncated normal, but I wouldn't have expected this behavior given that the posterior using the non-truncated normal wasn't ever dipping close to zero (or close to even the mean estimated by the truncated normal).

## Code

modfile<-nimbleCode( {

 mu.n ~ dunif(0.0001,8)

 sd~dinvgamma(0.001, 0.001)

  for (i in 1:T){

  data[i] ~ T(dnorm(mu.n,sd=sd),0,1)

#  data[i] ~ dnorm(mu.n,sd=sd)

 }

})

mod.inits <- list(mu.n = c(1), sd = c(0.5))

data<-c(0.013,0.077,0.069,0.055,0.176,0.03,0.097,0.097,0.108,0.006,0.016,0.016,0.125,0.015,0.004,0.013,0.022,0.042,0.032,0.035,0.017,0.031,0.032,0.564,0.01,0.023,0.017,0.123,0.387,0.097,0.147,0.11,0.265,0.114,0.003,0.053,0.062,0.012,0.087,0.034,0.066,0.075,0.162,0.017,0.121,0.276,0.01,0.575,0.205,0.008)

mod.data <- list(data = data)

mod.consts <- list(T = length(data))

results<-nimbleMCMC(code=modfile,constants=mod.consts,data=mod.data,inits=mod.inits,niter=50000,nburnin =25000)

plot(results[,1])

Thank you for any insights,

Franny

[Daniel Turek]

Franny, this is an interesting question.  I've looked at it for a while, and I believe all the inferences are correct, but it's only that the results appear non-intutitve for the truncated case.

For the non-truncated model, yes, the posterior distribution for mu.n nicely represents the mean of the data (both around 0.095, see code below).  That's what we expect.  But for the truncated model, the posterior of mu.n changes drastically, as you observed.  What is happening, is with the Normal likelihood truncated to the interval (0, 1), when the mean of this truncated likelihood is closer to 0 (in this case, centered around 0.0126), a large part of the of the Normal distribution falls below 0, and thus for the truncated likelihood, this larger fraction of probability density is re-distributed to the interval (0, 1), which serves to actually increase the overall likelihood of the data (as compared to having the mean of the truncated likelihood larger, as in the non-truncated case, which causes less of the Normal likelihood to be redistributed).  This effect is increased even further by the posterior distribution of the truncated Normal likelihood standard deviation (sd in your model) increasing slightly, to push more of the Normal likelihood probability below 0 (to be re-distributed in the interval (0, 1) and further increase the likelihood of the data).

You can see this in several ways.  The code below fits both models (truncated, non-truncated), then creates a plot to overlay two truncated likeloods (using the parameter posterior means fitted using each model) over a histogram of the data, so you can sort of "see" this effect of the truncated Normal likelihood.

To see this quantitatively, using the NIMBLE model object (modelTrunc) in the code below, where we look at the likelihood of the data, when setting the truncated Normal likelihood mean and standard deviation to the posterior mean parameter values from first the *non-truncated* model posteriors (giving a log-likelihood of the data of 46.2), and then using the posterior mean parameter values from the *truncated* model (giving a log-likelihood of the data of 55.6), verifying that in the truncated model, having the truncated Normal likelihood mean much closer to 0 does indeed increase the likelihood of the data.  I'll stress again, this result is because of the redistribution of the likelihood probability density, due to the truncation.

I hope this makes sense.  It was interesting to look into this, since I agree the main result seems non-intuitive.  The logical follow-up question is, how exactly you want to write your model given these observations?  Code is below.

Cheers,

Daniel

library(nimble)

## Code
modfile <- nimbleCode( {
    mu.n ~ dunif(0.0001,8)
    sd ~ dinvgamma(0.001, 0.001)

    for (i in 1:T){

        data[i] ~ dnorm(mu.n,sd=sd)
        ##data[i] ~ T(dnorm(mu.n,sd=sd),0,1)

    }
})

mod.inits <- list(mu.n = c(1), sd = c(0.5))

data <- c(0.013,0.077,0.069,0.055,0.176,0.03,0.097,0.097,0.108,0.006,0.016,0.016,0.125,0.015,0.004,0.013,0.022,0.042,0.032,0.035,0.017,0.031,0.032,0.564,0.01,0.023,0.017,0.123,0.387,0.097,0.147,0.11,0.265,0.114,0.003,0.053,0.062,0.012,0.087,0.034,0.066,0.075,0.162,0.017,0.121,0.276,0.01,0.575,0.205,0.008)

mod.data <- list(data = data)
mod.consts <- list(T = length(data))

set.seed(0)
results <- nimbleMCMC(code=modfile,constants=mod.consts,data=mod.data,inits=mod.inits,niter=50000,nburnin =25000)

modfileTrunc <- nimbleCode( {
    mu.n ~ dunif(0.0001,8)
    sd ~ dinvgamma(0.001, 0.001)

    for (i in 1:T){

        ##data[i] ~ dnorm(mu.n,sd=sd)
        data[i] ~ T(dnorm(mu.n,sd=sd),0,1)
    }
})

set.seed(0)
resultsTrunc <- nimbleMCMC(code=modfileTrunc,constants=mod.consts,data=mod.data,inits=mod.inits,niter=50000,nburnin =25000)

samplesSummary(results)
##            Mean     Median    St.Dev.  95%CI_low 95%CI_upp
## mu.n 0.09498151 0.09497028 0.01825259 0.05967935 0.1316005
## sd   0.12761111 0.12641033 0.01317482 0.10505308 0.1570089

mean(data)
## [1] 0.09502

samplesSummary(resultsTrunc)
##            Mean      Median    St.Dev.    95%CI_low  95%CI_upp
## mu.n 0.01261093 0.009554858 0.01114038 0.0004931849 0.04167774
## sd   0.15632670 0.154910845 0.01628674 0.1283629081 0.19204721


## plot histogram of the data, with truncated Normal likelihood
## curves, one using the non-truncated parameter posterior means, and the
## other using the truncated parameter posterior means
xlim <- c(0, 1)
ylim <- c(0, 10)
plot(-10, -10, xlim=xlim, ylim=ylim, xlab='', ylab='')
hist(data, freq = FALSE, breaks = 20, add = TRUE)
dnormTrun01 <- function(x, mu, s) dnorm(x, mu, s) / (pnorm(1, mu, s) - pnorm(0, mu, s))
plot(function(x) dnormTrun01(x, 0.01261093, 0.15632670), add = TRUE, col = 'blue')
plot(function(x) dnormTrun01(x, 0.09498151, 0.12761111), add = TRUE, col = 'red')
legend('topright', legend = c('Truncated', 'Non-truncated'), col = c('blue', 'red'), lwd=1)

modelTrunc <- nimbleModel(code=modfileTrunc,constants=mod.consts,data=mod.data,inits=mod.inits)

## use the non-truncated parameter posterior means:
modelTrunc$mu.n <- 0.09498151
modelTrunc$sd <- 0.12761111
modelTrunc$calculate()
## likelihood of the data:
modelTrunc$getLogProb('data')
## [1] 46.18509

## now, using the truncated parameter posterior means:
modelTrunc$mu.n <- 0.01261093
modelTrunc$sd <- 0.15632670
modelTrunc$calculate()
## likelihood of the data:
modelTrunc$getLogProb('data')
## [1] 55.60003    ## larger

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/cc979812-46e9-4135-8ba0-e520e2a53224n%40googlegroups.com.