new-sampling: specifying a new sampling distribution (illustrated with normal)

This example includes a demonstration of how future observations can also be modelled.

Model:

model {

   #                  trick using zeros

   for (i in 1:7) {

      zeros[i] <- 0

      zeros[i] ~ dpois(phi[i])         # likelihood is exp(-phi[i])

      phi[i] <- log(sigma[1]) + 0.5 * pow((x[i] - mu[1]) / sigma[1], 2)         # -log(likelihood)

   }

   #                  check using normal distribution

   for (i in 1:7) {

      x[i] ~ dnorm(mu[2], prec)

   }

   prec <- 1 / (sigma[2] * sigma[2])

   for (k in 1:2) {

      mu[k] ~ dunif(-10, 10)

      sigma[k] ~ dunif(0, 10)

   }

}

Data:

list(x = c(-1, -0.3, 0.1, 0.2, 0.7, 1.2, 1.7))

Initial values:

list(sigma = c(1, 1), mu = c(0, 0))

The results show a close match between the standard method and the zeros trick:

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample

   mu[1]   0.3661   0.4893   0.00715   -0.6296   0.364   1.336   1001   5000

   mu[2]   0.3668   0.4734   0.007003   -0.5591   0.3677   1.291   1001   5000

   sigma[1]   1.191   0.4902   0.01456   0.6293   1.075   2.476   1001   5000

   sigma[2]   1.168   0.461   0.009039   0.6265   1.062   2.311   1001   5000

Modelling future observations or missing data - add a missing data point to the data-file

Model:

model {

   #                  trick using zeros

   for (i in 1:8) {

      zeros[i] <- 0

      zeros[i] ~ dpois(phi[i])         # likelihood is exp(-phi[i])

      phi[i] <- log(sigma[1]) + 0.5 * pow((x[i] - mu[1]) / sigma[1], 2)         # -log(likelihood)

      x[i] ~ dflat()         # improper uniform prior on x's

   }

   #                  check using normal distribution

   for (i in 1:8) {

      x.rep[i] ~ dnorm(mu[2], prec)

   }

   prec <- 1 / (sigma[2] * sigma[2])

   for (k in 1:2) {

      mu[k] ~ dunif(-10, 10)

      sigma[k] ~ dunif(0, 10)

   }

}

Data:

list(x = c(-1, -0.3, 0.1, 0.2, 0.7, 1.2, 1.7, NA), x.rep=c(-1, -0.3, 0.1, 0.2, 0.7, 1.2, 1.7, NA))

(Note that in this case we need to replicate the dataset in order to compare results.)

Initial values:

list(sigma = c(1, 1), mu = c(0, 0), x = c(NA, NA, NA, NA, NA, NA, NA, 0))

(We need to give an initial value to the missing data, otherwise WinBUGS will try to generate one from the improper prior and crash!)

The agreement is fine but the MC error on the prediction is large and so a long run is necessary:

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample

   mu[1]   0.3689   0.4948   0.003145   -0.6157   0.3714   1.351   5001   50000

   mu[2]   0.368   0.486   0.002319   -0.6091   0.3713   1.326   5001   50000

   sigma[1]   1.194   0.5066   0.005429   0.6222   1.075   2.46   5001   50000

   sigma[2]   1.187   0.4989   0.003917   0.6224   1.069   2.45   5001   50000

   x[8]   0.357   1.407   0.01511   -2.494   0.3608   3.133   5001   50000

   x.rep[8]   0.3585   1.381   0.006126   -2.372   0.3584   3.105   5001   50000