Latent variable probit model

[Ethan]

I am trying to utilize the latent variable framework of Chib (1993) to exploit conjugacy for a probit regression. The basic idea of this latent variable representation is that if y ~ Bern(Phi(x’b)) then we can write

Pr(y = 1 | x, b) = \int Pr(y = 1 | y*, x, b) f(y* | x , b) dy*

where Pr(y = 1 | y*, x, b) = 1 if y* > 0 and 0 otherwise (i.e., y is deterministic given y*).

More importantly, y* | y = 1 ~ TN(x’b, 1, a = 0, b = Inf) and y* | y = 0 ~ TN(x’b, 1, a = -Inf, b = 0), which is easy to sample from.

The Nimble code I have is as follows:

code <- nimbleCode({
  for ( i in 1:n ) {
    mean[i]  <- beta0 + beta1 * x[i]
    ystar[i] ~ dnorm(mean[i], sd = 1)
    y[i]     ~ dinterval(ystar[i], 0)
  }
  beta0 ~ dnorm(0, sd = 10)
  beta1 ~ dnorm(0, sd = 10)
})

While this appears to work in terms of getting an answer close to the MLE, the sampler does not exploit the conjugacy of ystar (although, interestingly, it does use a conjugate sampler for beta).

I have attempted to write the quantile and random number generator function as follows:

qlatentnorm = nimbleFunction(
  run = function(
    p = double(0), mean = double(0, default = 0), y = integer(0)
    , lower.tail = integer(0, default = 1)
    , log.p = integer(0, default = 0)
  ) {
    returnType(double(0))
    if (log.p)
      p <- exp(p)
    if (!lower.tail)
      p <- 1 - p
    if ( y == 0 ) {
      logp <- log(p) + pnorm(0, mean = mean, log.p = 1)
      return( qnorm(logp, mean = mean, log.p = 1) )
    } else {
      lowProb <- pnorm(0, mean = mean)
      return( qnorm( p + lowProb * (1 - p), mean = mean ) )
    }
    return(0)  ## never reached
  }
)

rlatentnorm = nimbleFunction(
  run = function(
    n = integer(0), mean = double(0, default = 0), y = integer(0)
  ) {
    returnType(double(0))
    return( qlatentnorm(runif(1), mean = mean, y = y) )
  }
)

Both of these functions compile fine, and a histogram vs. density plot suggests that these functions have been implemented correctly.

I have programmed sampling the y*’s as follows:

latentNormalSampler <- nimbleFunction(
  name = 'latentNormalSampler',
  contains = sampler_BASE,
  setup = function(model, mvSaved, target, control) {
    ## Get mean parameters
    means <- model$getParents(target)
    ## Get corresponding binary variables
    ybin <- model$getDependencies(target, stochOnly = TRUE, self = FALSE)
    ## Get number of latent variables to update
    nupdates <- length(ybin)
  }
  , run = function() {
      ## Empty vector for draws
      draw <- rep(0, nupdates)
      for ( i in 1:nupdates ) {
        mean_i  <- values(model, means[i])[1]
        ybin_i  <- values(model, ybin[i])[1]
        draw[i] <- rlatentnorm(1, mean_i, ybin_i)
      }
      model[[target]] <<- draw
  },
  methods = list(
    reset = function() { }
  )
)

While this is lightning fast, it produces results that do not align at all with the MLE.

[Daniel Turek]

Ethan, thanks for the email, and great that you're getting into the programming aspects of nimble and MCMC samplers.

Without completely processing your sampling algorithm, can you try making one change (below) to your sampler, and seeing if the results are more what you'd expect?  Basically, after updating the value of the target node in the model, also (1) updating the associated log-density of [the new value of] the target node, using a call to model$calculate(target), and also (2; this is necessary for the MCMC to operate correctly on the entire model) updating the same values (the value of the target node, and the corresponding log-density) in the "mvSaved" data structure.  This is done using a call to nimCopy(...), and this agreement between the values in "model" and "mvSaved" objects is a pre-condition / post-condition requirement of all nimble MCMC samplers.  See code below, with 2 lines added.

latentNormalSampler <- nimbleFunction(
    name = 'latentNormalSampler',
    contains = sampler_BASE,
    setup = function(model, mvSaved, target, control) {
        ## Get mean parameters
        means <- model$getParents(target)
        ## Get corresponding binary variables
        ybin <- model$getDependencies(target, stochOnly = TRUE, self = FALSE)
        ## Get number of latent variables to update
        nupdates <- length(ybin)
    },
    run = function() {
        ## Empty vector for draws
        draw <- rep(0, nupdates)
        for ( i in 1:nupdates ) {
            mean_i  <- values(model, means[i])[1]
            ybin_i  <- values(model, ybin[i])[1]
            draw[i] <- rlatentnorm(1, mean_i, ybin_i)
        }
        model[[target]] <<- draw

        model$calculate(target)
        nimCopy(from = model, to = mvSaved, row = 1, nodes = target, logProb = TRUE)

    },
    methods = list(
        reset = function() { }
    )
)

All that said, I have not parsed the logic of your rlatentnorm function, nor the latentNormalSampler, but this did jump out at me.

Please let us know if this works.

Cheers,
Daniel

--
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/64742c0c-deb9-4e89-877b-64d582141a99n%40googlegroups.com.

[Ethan]

Thanks for the quick response, Daniel!

I have replaced my original code with yours, but it still does not give the expected results.

Regarding the functions, rlatentnorm exploits the fact that if X ~ F, then we can generate from F via F^(-1)(U) where U ~ Unif(0, 1). Thus, provided the quantile function is correct (which I checked), the rlatentnorm function should work as expected.

Regarding latentNormalSampler, what I am trying to do is replace the default RW_sampler for ystar with the conjugacy result. If y == 0, then ystar ~ TruncatedNormal(x'beta, 1, a = -Inf, b = 0) and if y == 1 then ystar ~ TruncatedNormal(x'beta, 1, a = 0, b = Inf)

Honestly, it would be easier to probably just say:

draw[i] ~ (2 * ybin[i] - 1) * abs(rnorm(1))

But I don't think the issue is with the random number generation function.

For completeness, here's a minimally reproducible example

library(nimble)
library(dplyr)
library(posterior)
library(bayesplot)

dlatentnorm <- nimbleFunction(
  run = function(
    x = double(0), mean = double(0, default = 0), y = integer(0), log = integer(0, default = 0)
  ) {
    returnType(double(0))
    if ( x < 0 & y == 1 ) {
      logProb <- -Inf
    } else if ( x > 0 & y == 0) {
      logProb <- -Inf
    } else {
      logProb <- dnorm(x, mean = mean, sd = 1, log = 1) - pnorm(0, mean = mean, sd = 1, lower.tail = (1-y), log.p = 1)
    }
    if ( log )
      return(logProb)
    return(exp(logProb))
  }
)
platentnorm <- nimbleFunction(
  run = function(
    q = double(0), mean = double(0, default = 0), y = integer(0)

    , lower.tail = integer(0, default = 1)
    , log.p = integer(0, default = 0)
  ) {
    returnType(double(0))

    if ( y == 0 ) {

      if ( q >= 0 ) {
        logProb <- Inf
      } else {
        logProb <- pnorm(q, mean, 1, log.p = 1) - pnorm(0, mean, 1, log.p = 1)
      }
    } else {
      if ( q < 0 ) {
        logProb <- -Inf
      } else {
        lowerProb <- pnorm(0, mean = mean, sd = 1)
        logProb <- log( pnorm(q, mean) - lowerProb ) - log(1 - lowerProb)
      }
    }
    if ( lower.tail ) {
      if ( log.p )
        return(logProb)
      else
        return(exp(logProb))
    } else {
      prob <- 1 - exp(logProb)
      if ( log.p ) {
        return(log(prob))
      } else {
        return(prob)

cdlatentnorm <- compileNimble(dlatentnorm)
cplatentnorm <- compileNimble(platentnorm)
cqlatentnorm <- compileNimble(qlatentnorm)
crlatentnorm <- compileNimble(rlatentnorm)

code <- nimbleCode({
  for ( i in 1:n ) {
    mean[i]  <- beta0 + beta1 * x[i]
    ystar[i] ~ dnorm(mean[i], sd = 1)
    y[i]     ~ dinterval(ystar[i], 0)
  }
  beta0 ~ dnorm(0, sd = 10)
  beta1 ~ dnorm(0, sd = 10)
})

## Simulate data
n     = 100
x     = rnorm(n, sd = 2)
ystar = rnorm(n, -1 + x, sd = 1)
y     = ifelse(ystar < 0, 0, 1)

## Run nimble code
const  <- list('n' = n)
data   <- list('y' = y, 'x' = x)
init   <- list('beta0' = 0, 'beta1' = 0, 'ystar' = rnorm(const$n, ifelse(data$y == 0, -1, 1), sd = 0.2) )
model  <- nimbleModel(code, constants = const, data = data, inits = init)
cmodel <- compileNimble(model)
conf   <- configureMCMC(model, monitors = c('beta0', 'beta1'))
conf$removeSampler('ystar')
conf$addSampler('ystar', 'latentNormalSampler')
mcmc   <- buildMCMC(conf)
cmcmc  <- compileNimble(mcmc, project = model)
smpl   <- runMCMC(cmcmc)
summ   <- smpl %>% summarize_draws()
summ
summary(glm(y ~ x, family = binomial('probit')))

[Daniel Turek]

Ethan, thanks for sending the reproducible example.  I've looked over your code, and I don't spot any error offhand.  I'll continue trying to look over it, but in the meantime, I suggest you try running your function uncompiled using runMCMC(mcmc), and using R's debugger tools.  This would let you step though / into the functions, and make certain they're operating as you expect.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/7a27de0c-aa2c-4c1f-b6b8-59639ce0a631n%40googlegroups.com.

[Chris Paciorek]

Hi Ethan,

it looks like the problem here is that the result of getParents is not what you are expecting -- it includes the `beta`s as well as the `mean`s:

Browse[2]> means
  [1] "beta0"     "beta1"     "mean[1]"   "mean[2]"   "mean[3]"   "mean[4]"  
  [7] "mean[5]"   "mean[6]"   "mean[7]"   "mean[8]"   "mean[9]"   "mean[10]"
 [13] "mean[11]"  "mean[12]"  "mean[13]"  "mean[14]"  "mean[15]"  "mean[16]" 

Take a look at the help information on getParents: `determOnly = TRUE` should get you what you need.

-chris

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CAKbe0hrVxEC%2BDTfJD1%3Dt9JF286bhtC%2BKQPTavTK6jq%3D%3D4L3MgQ%40mail.gmail.com.

[Ethan]

That seems to have worked. Thanks, Chris!