How to create a conjugate sampler for binomial-truncated beta

[Ethan]

Hello,

I need to run a model consisting of Bernoulli latent variables. I want to elicit a truncated prior on the latent variables, which is conjugate. However, if I try

gamma ~ T(dbeta(2, 2), 0, 0.5)

the resulting sampler does not recognize conjugacy.

I have written out the density, CDF, quantile function, and random number generation for a truncated beta, but I'm not sure how to create the sampler.

I've made some sample code to illustrate that it does not recognize conjugacy for a simple model (see below). 

## If Y = X 1{a < X < b}. Let f = PDF of Y and g/G = PDF/CDF of X.
## f(y) = 1 / [G(b) - G(a)] g(y) 1{a < y < b}
dtruncbeta <- nimbleFunction(
  run = function(
    x = double(0), shape1 = double(0), shape2 = double(0),
    a = double(0,default = 0), b = double(0, default = 1),
    log = integer(0, default = 0)
  ) {
    returnType(double(0))
    if ( x <= a ) {
      lp <- -Inf
    } else if ( x >= b ) {
      lp <- -Inf
    } else {
      lp      <- dbeta(x, shape1, shape2, log = TRUE)   ## initialize to beta PDF
      logCDFa <- -Inf
      logCDFb <- 0
      if ( a > 0 ) {
        logCDFa <- pbeta(a, shape1, shape2, log.p = TRUE)
      }
      if ( b < 1 ) {
        logCDFb <- pbeta(b, shape1, shape2, log.p = TRUE)
      }
      lognc <- logCDFb + log( 1 - exp(logCDFa - logCDFb) )
      lp    <- -lognc + lp
    }
    if(log){
      return(lp)
    }
    return(exp(lp))
  }
)


## If Y = X 1{a < X < b}. Let F = CDF of Y and G = CDF of X.
## F(y) = [G(y) - G(a)] / [G(b) - G(a)]
ptruncbeta <- nimbleFunction(
  run = function(
    q = double(0), shape1 = double(0), shape2 = double(0),
    a = double(0,default = 0), b = double(0, default = 1),
    log.p = integer(0, default = 0)
  ) {
    returnType(double(0))
    if ( x <= a ) {
      lp <- -Inf
    } else if (x >= b) {
      lp <- 0
    } else {
      logCDFx <- pbeta(x, shape1 = shape1, shape2 = shape2, log.p = TRUE)
      logCDFa <- -Inf
      logCDFb <- 1
      if ( a > 0 ){
        logCDFa <- pbeta(a, shape1, shape2, log.p = TRUE)
      }
      if ( b < 1 ) {
        logCDFb <- pbeta(b, shape1, shape2, log.p = TRUE)
      }
      logNum <- logCDFx + log(1 - exp(logCDFa - logCDFx))
      logDen <- logCDFb + log(1 - exp(logCDFa - logCDFb))
      lp     <- logNum - logDen
    }
    if(log.p) {
      return(lp)
    }
    return(exp(lp))
  }
)


qtruncbeta <- nimbleFunction(
  run = function(
    p = double(0), shape1 = double(0), shape2 = double(0),
    a = double(0,default = 0), b = double(0, default = 1),
    log.p = integer(0, default = 0)
  ) {
    returnType(double(0))
    if(log.p) {
      x <- exp(x)
    }
    if ( x <= 0 ) {
      return(NaN)
    } else if ( x >= 1) {
      return(NaN)
    } else {
      CDFa  <- pbeta(a, shape1, shape2, log.p = FALSE)
      CDFb  <- pbeta(b, shape1, shape2, log.p = FALSE)
      return( qbeta( CDFa + (CDFb - CDFa) * x , shape1 = shape1, shape2 = shape2) )
    }
    return(NaN)
  }
)

## Use inverse-CDF method for random number generation
rtruncbeta <- nimbleFunction(
  run = function(
    n = integer(0), shape1 = double(0), shape2 = double(0),
    a = double(0,default = 0), b = double(0, default = 1)
  )
    return(qtruncbeta(runif(n), shape1 = shape1, shape2 = shape2, a = a, b = b, log.p = FALSE))
)


n = 100
code <- nimbleCode({
  gamma ~ dtruncbeta(2, 2, 0, 0.75)
  # gamma ~ T(dbeta(2, 2), 0, 0.75)
  for ( i in 1:n ) {
    y[i] ~ dbinom(size = 1, prob = gamma)
  }
})

## Generate binomial data
y <- rbinom(100, size = 1, prob = 0.25)

## Model / configure
mod <- nimbleModel(code, constants = list('n' = n), data = list('y' = y), init = list('gamma' = 0.1))
modConf <- configureMCMC(mod)

[Daniel Turek]

Ethan, this is nice work.  Way to write the rtruncbeta, dtruncbeta, and p and q functions, for starters.  And moreso, way to write the density evaluations in a numerically stable way.  That's something we often run into problems with, and you seem to have addressed it well.

(I feel like there was a user discussion about adding / detecting conjugacy fairly recently, but on a quick search I just can't find it, so we get to start afresh).

Quickly noting also, I tried your code, and went through registerDistributions to register your distribution with Nimble (good measure):

registerDistributions(
    list(dtruncbeta = list(
             BUGSdist = 'dtruncbeta(shape1, shape2, a, b)',
             types = c('shape1 = double()', 'shape2 = double()', 'a = double()', 'b = double()'),
             pqAvail = TRUE)
         )
)

And it identifies the problem that the qtruncbeta and ptruncbeta need the final argument of lower.tail, which needs to be added as:

## If Y = X 1{a < X < b}. Let F = CDF of Y and G = CDF of X.
## F(y) = [G(y) - G(a)] / [G(b) - G(a)]
ptruncbeta <- nimbleFunction(
    run = function(
                   q = double(0), shape1 = double(0), shape2 = double(0),
                   a = double(0,default = 0), b = double(0, default = 1),

                   lower.tail = double(0, default = 0),

                   log.p = integer(0, default = 0)
                   ) {

qtruncbeta <- nimbleFunction(
    run = function(
                   p = double(0), shape1 = double(0), shape2 = double(0),
                   a = double(0,default = 0), b = double(0, default = 1),

                   lower.tail = double(0, default = 0),

                   log.p = integer(0, default = 0)
                   ) {

And then of course also using the lower-tail argument correctly in the implementation of both these functions.

Moving on, Nimble has a predefined list of conjugacies that it "detects".  You can look at them - which are sorted by the prior distribution, and for each prior, the "detectable conjugate likelihood distributions" follow, along with a semi-clear-semi-cryptic definition of the resulting conjugate posterior distribution (from which to draw samples).  Please take a quick look at this predefined list the "conjugacyRelationshipsInputList" at the very top of MCMC_conjugacy.R:

https://github.com/nimble-dev/nimble/blob/devel/packages/nimble/R/MCMC_conjugacy.R

Now, for nimble to automatically "detect" your conjugacy, this list object would need an element defining the conjugacy relationship for "dtruncbeta" with the dbinom prior.  This new list element would need to be written to correctly "define" the conjugate posterior, similar to in the other list entries.  I could talk you through this "semi cryptic way of defining the posterior" if that's helpful.

I see three options here:

(1) This would work, without too much trouble.  Create your own branch of the nimble repository.  On that branch, add a new element to the "conjugacyRelationshipsInputList", I did a quick draft of what the new element would look like, I'm not saying the conjugate posterior is correct, but the "scaffolding" is there for you (and it might be correct, just the same as beta-binomial conjugacy, but you should confirm that):

conjugacyRelationshipsInputList <- list(

    ## dtruncbeta
    list(prior = 'dtruncbeta',
         link = 'identity',
         dependents = list(
             dbinom = list(param = 'prob', contribution_shape1 = 'value', contribution_shape2 = '1 - value')),
         posterior = 'dtruncbeta(shape1 = prior_shape1 + contribution_shape1,
                                 shape2 = prior_shape2 + contribution_shape2,
                                 a      = prior_a,
                                 b      = prior_b)'),

ETC.... rest of the list

Then build the Nimble package using your modified branch, and it should work out-of-the-box.

(2) Forget about having Nimble "detect" it automatically for you.  Instead, just write your own conjugate sampler.  This would require writing a nimbleFunction, and getting a little into the Nimble DSL, but there's an example in the User Manual (section 15.5), and that could get you started.  Then, once you've written your sampler, as:

sampler_conjugate_dtrunbeta <- nimbleFunction(......)

You can easily add it to the MCMC configuration:

conf$removeSamplers('gamma')

conf$addSampler(target = 'gamma', type = 'sampler_conjugate_dtrunbeta')

And the rest will work.

(3) We, the Nimble development team, could consider extending the scope of where and how Nimble looks for conjugate relationships, so that there would be some function, say "registerConjugacy", which allows users to register their own (new) conjugate relationship.  This would be similar to registerDistributions.  This is something we should probably consider adding, it sounds like, but it's not going to happen today.  We'll discuss it, though.

Anyway, I hope this gets you going.  Keep in touch.

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/9d516e47-fca2-42c7-8f33-4683d91439ddn%40googlegroups.com.

[Ethan]

Hi Daniel,

Thanks a lot! I would like to try just writing my own sampler, but I think I need a little more than what is provided in Section 15.5. For the RW Metropolis example, one doesn't really need access to the other parameters / data than the one being considered.

How could I access the number of binomial success and failures?

[Daniel Turek]

Ethan, good question.  The code below (adapted from your original code) should get you started no problem.  Please make sure to take a look at the comments.

Noting that this "mySampler" doesn't actually do anything, but maybe doing anything is overrated anyway.

library(nimble)

n = 10

code <- nimbleCode({

    gamma ~ T(dbeta(2, 2), 0, 0.75)
    for ( i in 1:n ) {
        y[i] ~ dbinom(size = 1, prob = gamma)
    }
})

set.seed(0)

## Generate binomial data
y <- rbinom(n, size = 1, prob = 0.25)

## Model / configure
mod <- nimbleModel(code, constants = list('n' = n), data = list('y' = y), init = list('gamma' = 0.1))

mySampler <- nimbleFunction(
    name = 'mySampler',
    contains = sampler_BASE,
    setup = function(model, mvSaved, target, control) {
        ## query the model, to get all the stochastic dependents
        ## of the 'target' node:
        dependentNodes <- model$getDependencies(target, stochOnly = TRUE, self = FALSE)
        ## just for the for-loop in the run code:
        nDependentNodes <- length(dependentNodes)
        ## do a quick check here, that all
        ## the stochastic dependents are 'dbin'-distributed.
        ## note NIMBLE has changed the distribution name to 'dbin'.
        ## but you could put whatever checks you want here, or none.
        if(!all(model$getDistribution(dependentNodes) == 'dbin')) stop('dependent nodes are not binomial')
    },
    run = function() {
        for(i in 1:nDependentNodes) {

            ## get parameters of the dependent node distributions:
            size <- model$getParam(dependentNodes[i], 'size')
            prob <- model$getParam(dependentNodes[i], 'prob')

            ## get the value, number of successes or failures, of each dependent node:
            value <- values(model, dependentNodes[i])
            print(dependentNodes[i], ':')
            print('size=', size, ', prob=', prob, ', value=', value)
        }
        ## Note that you can also use:
        ## values(model, dependentNodes)
        ## to get a vector containing *all* the dependent node values:
        valuesVector <- values(model, dependentNodes)
        print('all dependent node values:')
        print(valuesVector)
    },
    methods = list(
        reset = function() { }
    )
)

modConf <- configureMCMC(mod)

modConf$removeSamplers('gamma')

modConf$addSampler(target = 'gamma', type = 'mySampler')

## print sampler configuration:
modConf

## build MCMC:
Rmcmc <- buildMCMC(modConf)

## compile model and MCMC:
Cmod <- compileNimble(mod)
Cmcmc <- compileNimble(Rmcmc, project = mod)

## run:
Rmcmc$run(1)   ## uncompiled
Cmcmc$run(1)   ## compiled

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/b8c0dcf8-905a-4eef-a603-bad2388f9442n%40googlegroups.com.

[Ethan]

Great! I was able to figure out how to complete the sampler. As a small note, the values function, when compiled, returns a vector and threw an error. I thus replaced:   

value <- values(model, dependentNodes[i])

with

value <- values(model, dependentNodes[i])[1]

which seems to work (as an aside, how do you get fixed width text in Google groups?)

The following code appears to work, checking against the known conjugacy result (in the univariate case). Feel free to post on as an example if you wish. I appreciate the help!

set.seed(123)

## If Y = X 1{a < X < b}. Let f = PDF of Y and g/G = PDF/CDF of X.
## f(y) = 1 / [G(b) - G(a)] g(y) 1{a < y < b}
dtruncbeta <- nimbleFunction(
  run = function(
    x = double(0), shape1 = double(0), shape2 = double(0),

    a = double(0, default = 0), b = double(0, default = 1),

      p <- exp(p)
    }
    if ( p <= 0 ) {
      return(NaN)
    } else if ( p >= 1) {

      return(NaN)
    } else {
      CDFa  <- pbeta(a, shape1, shape2, log.p = FALSE)
      CDFb  <- pbeta(b, shape1, shape2, log.p = FALSE)
      return(

        qbeta( CDFa + (CDFb - CDFa) * p , shape1 = shape1, shape2 = shape2)
      )
    }
  }

)

## Use inverse-CDF method for random number generation
rtruncbeta <- nimbleFunction(
  run = function(
    n = integer(0), shape1 = double(0), shape2 = double(0),
    a = double(0,default = 0), b = double(0, default = 1)

  ){
    returnType(double(0))
    return(qtruncbeta(runif(1), shape1 = shape1, shape2 = shape2, a = a, b = b, log.p = FALSE))
  }
)

## Generate binomial data
n = 100

y <- rbinom(n, size = 1, prob = 0.25)

## Custom MCMC sampler for truncated beta
conjTruncBeta <- nimbleFunction(
  name = 'conjTruncBeta',

  contains = sampler_BASE,
  setup = function(model, mvSaved, target, control) {
    ## query the model, to get all the stochastic dependents
    ## of the 'target' node:
    dependentNodes <- model$getDependencies(target, stochOnly = TRUE, self = FALSE)
    ## just for the for-loop in the run code:
    nDependentNodes <- length(dependentNodes)

  },
  run = function() {
    shape1_post <- model$getParam(target, 'shape1')    ## initialize to prior
    shape2_post <- model$getParam(target, 'shape2')    ## initialize to prior

    for(i in 1:nDependentNodes) {
      ## get parameters of the dependent node distributions:

      size        <- model$getParam(dependentNodes[i], 'size')   ## size of y[i]
      value       <- values(model, dependentNodes[i])[1]         ## num successes in y[i]
     
      ## Compute posterior quantities based on conjugacy
      shape1_post <- shape1_post + value
      shape2_post <- shape2_post + (size - value)
    }
    ## Draw from truncated beta using conjugacy
    draw  <- rtruncbeta(
      1, shape1 = shape1_post, shape2 = shape2_post,
      a = model$getParam(target, 'a'),
      b = model$getParam(target, 'b')
    )
    ## Store draw in model
    model[[target]] <- draw
   
    ## Store draw in mvsaved
    copy(from = model, to = mvSaved, row = 1,
         nodes = dependentNodes, logProb = TRUE)

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

## Code for model
code <- nimbleCode({
  gamma ~ dtruncbeta(2, 2, 0.10, 0.75)

  for ( i in 1:n ) {
    y[i] ~ dbinom(size = 1, prob = gamma)
  }
})

## Configure the model and MCMC and compile both
mod  <- nimbleModel(

  code,
  constants = list('n' = n),
  data = list('y' = y),

  init = list('gamma' = 0.2)
)
mod.comp <- compileNimble(mod)
mod.conf <- configureMCMC(mod)
mod.conf$removeSampler('gamma')
mod.conf$addSampler('gamma', 'conjTruncBeta')
mod.mcmc      <- buildMCMC(mod.conf)
mod.mcmc.comp <- compileNimble(mod.mcmc, project = mod)
smpl          <- runMCMC(mod.mcmc.comp, niter = 100000)

## Compare with direct sampling
nsucc <- sum(y)
nfail <- n - nsucc
hist(smpl, freq = F)
x  <- seq(0.10, 0.75, length.out = 5000)
fx <- sapply(x, function(x) dtruncbeta(x, nsucc + 2.0, nfail + 2.0, a = 0.10, b = 0.75))
lines(x = x, y = fx, col = 'blue')