BUGS Implementation Tricks

Implementing Custom Distributions Without Low-Level Coding in Previous Version of BUGS

In JuliaBUGS, users can simply create new distributions using the Distributions.jl interface and use them as built-in distributions. In previous versions of BUGS, defining new distributions required lower-level programming, as users needed to work directly with the underlying implementation.

Here we present some of the tricks that were used in previous BUGS implementations to create custom distributions without implementing them directly at the low level. These approaches are still valid in JuliaBUGS, though the native distribution interface is generally preferred.

The "Zeros Trick"

When you need a distribution not included in the standard set, the "zeros trick" offers an elegant solution. This technique leverages the fact that a Poisson observation with mean $\phi$ and value 0 has likelihood $e^{-\phi}$.

How it works:

Create artificial data points of zeros
Set $\phi[i] = -\log(L[i]) + C$, where:
- \[L[i]\]
  is your desired likelihood term
- \[C\]
  is a constant large enough to ensure $\phi[i] > 0$

Implementation:

C <- 10000    # Large constant ensuring phi[i] > 0

for (i in 1:N) {
    zeros[i] <- 0
    phi[i] <- -log(L[i]) + C
    zeros[i] ~ dpois(phi[i])
}

This method is particularly useful for implementing truncated distributions or any arbitrary likelihood function.

The "Ones Trick"

An alternative approach uses Bernoulli observations fixed at 1:

How it works:

Create artificial data points of ones
Define probabilities proportional to your desired likelihood: $p[i] = \frac{L[i]}{C}$
Choose $C$ large enough to ensure $p[i] < 1$

Implementation:

C <- 10000    # Large constant ensuring p[i] < 1

for (i in 1:N) {
    ones[i] <- 1
    p[i] <- L[i] / C
    ones[i] ~ dbern(p[i])
}

Using `dloglik` in OpenBUGS and MultiBUGS

The dloglik distribution provides a more direct approach for implementing custom likelihoods:

dummy[i] <- 0
dummy[i] ~ dloglik(logLike[i])

Where logLike[i] is the log-likelihood contribution for observation $i$. This essentially implements the "zeros trick" behind the scenes.

Example: Manual Normal Likelihood Implementation

model {
   for (i in 1:7) {
      dummy[i] <- 0
      dummy[i] ~ dloglik(logLike[i])
      logLike[i] <- -log(sigma) - 0.5 * pow((x[i] - mu)/sigma, 2)         
   }
   mu ~ dunif(-10, 10)
   sigma ~ dunif(0, 10)
}

Standard equivalent:

model {
   for (i in 1:7) {
      x[i] ~ dnorm(mu, prec)
   }
   prec <- 1 / (sigma * sigma)
   mu ~ dunif(-10, 10)
   sigma ~ dunif(0, 10)
}

Implementing Custom Prior Distributions

You can use dloglik to implement non-standard prior distributions:

theta ~ dflat()           # Use flat improper prior as base
dummy <- 0
dummy ~ dloglik(logLike)  # Add custom prior via log-likelihood
logLike <- log(desired_prior_for_theta)

Example: Manual Normal Prior Implementation

model {
   for (i in 1:7) {
      x[i] ~ dnorm(mu, prec)
   }
   dummy <- 0
   dummy ~ dloglik(phi)
   phi <- -0.5 * pow(mu, 2)  # log(N(0,1))
   mu ~ dflat()              # Base distribution
   prec <- 1 / (sigma * sigma)
   sigma ~ dunif(0, 10)
}

Standard equivalent:

model {
   for (i in 1:7) {
      x[i] ~ dnorm(mu, prec)
   }
   mu ~ dnorm(0, 1)
   prec <- 1 / (sigma * sigma)
   sigma ~ dunif(0, 10)
}

Note: Using dloglik for priors may trigger Metropolis sampling, potentially leading to slower convergence and higher Monte Carlo errors.

Working with Predictions and Complex Models

Predicting New Observations

To generate predictions for a new observation x.pred, specify it as missing and assign an improper uniform prior:

x.pred ~ dflat()  # Improper uniform prior

Be aware this approach may increase computational inefficiency and Monte Carlo error.

Handling Model Mixtures of Different Complexity

For mixture models with components of varying complexity, a standard mixture distribution approach is often sufficient without requiring reversible jump techniques:

model {
   mu ~ dunif(-5, 5)
   p ~ dunif(0, 1)
   m[1] <- 0       # First component mean
   m[2] <- mu      # Second component mean
   
   for (i in 1:100) {
      group[i] ~ dbern(p)           # Component membership
      index[i] <- group[i] + 1
      y[i] ~ dnorm(m[index[i]], 1)  # Observation from selected component
   }
}

Managing Random Set Sizes

When loop bounds depend on random quantities (e.g., changepoints), use step functions to conditionally include elements:

for (i in 1:N) {
   ind[i] <- 1 + step(i - K - 0.01)  # 1 if i ≤ K, 2 if i > K
   y[i] ~ model[ind[i]]              # Select appropriate model
}

Example: Computing the Sum of First K Integers

model {
   # Define possible values for K
   for (i in 1:10) {
      p[i] <- 1/10  # Equal probability for each value
      x[i] <- i     # Value i
   }
   
   # Random selection of K
   K ~ dcat(p[])
   
   # Sum elements conditionally
   for (i in 1:10) {
      xtosum[i] <- x[i] * step(K - i + 0.01)  # Include x[i] only if i ≤ K
   }
   
   # Compute final sum
   s <- sum(xtosum[])
}