Mixture Model using R2WinBUGS

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Merlise Clyde.Lecture Notes Mixture Model using WinBUGS.pdf (171.01 KB, 需要: 5 个论坛币)

2013-12-7 06:31:22 上传

Mixture Model using R2WinBUGS
需要: 5 个论坛币  [购买]

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：R2WinBUGS winbugs Mixture WINBUG Using

################
#  WinBUGS code for Finite Mixture model (Carlin and Louis, 2008, Example 2.19)
#  data are "eyes" data from WinBUGS manual
################

model
{
for( i in 1 : N ) {
   y[i] ~ dnorm(mu[i], tau)
   mu[i] <- lambda[T[i]]
   T[i] ~ dcat(P[])
   }
P[1:2] ~ ddirch(alpha[])
theta ~ dnorm(0.0, 1.0E-6) I(0.0, )
lambda[2] <- lambda[1] + theta
lambda[1] ~ dnorm(0.0, 1.0E-6)

Ystar1 ~ dnorm(lambda[1], tau)
Ystar2 ~ dnorm(lambda[2], tau)
Tstar ~ dcat(P[]);  astar <- Tstar-1
Ystar <- Ystar1*(1-astar) + Ystar2*astar

sigma ~ dunif(0.01,100)   #  vague Gelman prior for sigma
tau<- 1/(sigma*sigma)
}

DATA:
list(y = c(529.0, 530.0, 532.0, 533.1, 533.4,
           533.6, 533.7, 534.1, 534.8, 535.3,
           535.4, 535.9, 536.1, 536.3, 536.4,
           536.6, 537.0, 537.4, 537.5, 538.3,
           538.5, 538.6, 539.4, 539.6, 540.4,
           540.8, 542.0, 542.8, 543.0, 543.5,
           543.8, 543.9, 545.3, 546.2, 548.8,
           548.7, 548.9, 549.0, 549.4, 549.9,
           550.6, 551.2, 551.4, 551.5, 551.6,
           552.8, 552.9,553.2), N = 48, alpha = c(1, 1),
    T = c(1, NA, NA, NA, NA, NA, NA, NA, NA, NA,
          NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
          NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
          NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
          NA, NA, NA, NA, NA, NA, NA, 2))


INITS (two sets):

list(lambda = c(535, NA), theta = 5, sigma = 1.0, P=c(0.5, 0.5))
list(lambda = c(535, NA), theta =15, sigma = 5.0, P=c(0.2, 0.8))

Shifted-Wald mixture models in WinBUGS


Suppose you want a two component normal model. You might use something like this:

    model {
      for(i in 1:nobs){
        C ~ dcat(pi[1:2])
        d ~ dnorm(mu[C],tau[C])
      }
      pi[1:2] ~ ddirch(xi[1:2])
      for(j in 1:2){
        xi[j] <- 1
      }
      mu[1] ~ dnorm(1,1)
      tau[1] ~ dgamma(1,.1)
      mu[2] ~ dnorm(5,1)
      tau[2] ~ dgamma(4,.1)
      }
    }
In this model, C indexes the mixture components, and is distributed as a categorical random variable, which is to say that it takes the value ’1′ with probability pi[1] and the value ’2′ with probability pi[2] = 1-pi[1]. The vector pi is distributed as a Dirichlet random variable with parameters 1 and 1, which is to say that pi[1] is uniformly distributed. The vector d has the data in it, and the idea is that each element of d is being modeled as a normal variate with either mu[1] and tau[1] or mu[2] and tau[2], each of which is given its own distinct prior.

This model works just fine, so you might think that, after getting WinBUGS to recognize the shifted Wald distribution, it would work to make the appropriate changes and, say, try this (the parameters in the priors are arbitrary here – I’m just trying to get the model to work at all right now):

    model {
      for(i in 1:nobs){
        C ~ dcat(pi[1:2])
        d ~ ShiftedWald(v[C],a[C],b[C])
      }
      pi[1:2] ~ ddirch(xi[1:2])
      for(j in 1:2){
        xi[j] <- 1
      }
      v[1] ~ dgamma(8,4)
      a[1] ~ dgamma(10,2)
      b[1] ~ dgamma(6,3)
      v[2] ~ dgamma(2,1)
      a[2] ~ dgamma(4,4)
      b[2] ~ dgamma(16,2)
      }
    }
But this doesn’t work, and I have no idea why. When I try to run this, I get an ‘index out of range’ trap error. To the (small) extent that I can figure out what’s going wrong, the dcat distribution is allowing C to take values other than 1 or 2. Now, because I only want two mixture components, I can do the following instead, which seems to be working just fine (on simulated data):

    model {
      for(i in 1:nobs){
        C ~ dbern(pi)
        Cp <- C + 1
        d ~ ShiftedWald(v[Cp],a[Cp],b[Cp])
      }
      pi ~ dbeta(xi,nu)
      xi <- 1
      nu <- 1
      v[1] ~ dgamma(8,4)
      a[1] ~ dgamma(10,2)
      b[1] ~ dgamma(6,3)
      v[2] ~ dgamma(2,1)
      a[2] ~ dgamma(4,4)
      b[2] ~ dgamma(16,2)
      }
    }

Eyes: normal mixture models

http://www.mrc-bsu.cam.ac.uk/bug ... xampVol2/node5.html

Hierarchical Spatial N-mixture and Occupancy Models using WinBUGS

http://www.esapubs.org/archive/ecol/E092/036/suppl-1.htm

Fitting N-mixture models with random observer effects in BUGS and ADMB

BayesMix: An R package for Bayesian Mixture Modeling

Bettina Grun
Department of Statistics and Probability Theory
Vienna University of Technology

Abstract
The R package BayesMix provides the functionality for estimating univariate Gaussian mixture models with MCMC methods. Within a given model class users can modify the prior specifications and the initial values for developing a suitable model for their data. Furthermore, tools for analyzing the output of the MCMC simulations as, e.g., diagnostic plots, are available. The package is intended to be accompanying material for Fruhwirth-Schnatter(2006). With this package the user can comfortably reproduce some of the results presented in the book. The data sets and functions for generating the initial values and prior

specifications in the book are provided.

http://statmath.wu.ac.at/~gruen/BayesMix/bayesmix-intro.pdf

Logistic Regression using R2WinBUGS
http://d.hatena.ne.jp/trashmosh/20120623/p3

Bayesian Generalized Linear Models using R2WinBUGS

https://aquila2.iseg.utl.pt/aqui ... tFile&fileId=132907

model
{
for( i in 1 : N ) {
r ~ dbin(p,n)
logit(p) <- alpha + beta * (x - mean(x[]))
rhat <- n * p
phat <- r/n
}
beta ~ dnorm(0.0, 0.001)
alpha ~ dnorm(0.0, 0.001)
}
> data.b <- c("x", "n", "r", "N"); par.b <- c("alpha", "beta")
> bin.reg <- bugs(data.b, inits=NULL, par.b, "model-beetles1.txt",
n.chains = 3, n.iter = 20000, n.thin=1, bugs.directory = bugsdir,
+ working.directory = getwd(), clearWD=TRUE, debug=TRUE)
>
> print(bin.reg, digits=3)

alpha <- bin.reg$sims.array[,1 ,"alpha"]
beta <- bin.reg$sims.array[,1 ,"beta"]
invlogit <- function(x)
{
1/(1 + exp(-x))
}
plot(x, r/n, cex=2, col="red", pch="o", type="b")
for(b in 1: 50){
points(x, invlogit((alpha - beta*mean(x)) + beta * x),
cex=0.5, type="l", col="grey", lty=3)
}

Simple linear regression using R2OpenBUGS

Simple linear regression using R2OpenBUGS.pdf (11.33 MB, 需要: 5 个论坛币)

2013-12-7 11:42:07 上传

Simple linear regression using R2OpenBUGS
需要: 5 个论坛币  [购买]

First, let's load the package:
library(R2OpenBUGS)Next, write the model as though it was an R function. Note that the function does not have any arguments. Syntax used here is BUGS syntax. The function will be translated into BUGS code, so it's important to retain appropriate syntax for BUGS. If using constraints on parameters e.g. “I(0,)” then the syntax differs from the standard BUGS syntax - see the later example (Alternative 4). Note that in BUGS the ordering of likelihood and priors does not matter. The fact that we use alpha and beta before defining the prior does not matter. The convention is to specify the likelihood first and priors second.

linemodel <- function() {   
for (j in 1:N) {        Y[j ~ dnorm(mu[j, tau)  ## Response values Y are Normally distributed        
mu[j <- alpha + beta * (x[j - xbar)  ## linear model with x values centred   
}   

## Priors   
alpha ~ dnorm(0, 0.001)   
beta ~ dnorm(0, 0.001)   
tau ~ dgamma(0.001, 0.001)   
sigma <- 1/sqrt(tau)}

Next, let's prepare the data for analysis. BUGS accepts named list format data, including data structures with named columns (or mixtures of these). Note that BUGS doesn't react well if data passed to it is not used in the model… The best advice is to only pass what you need. Fortunately R is the perfect tool for preparing data for BUGS.

linedata <- list(Y = c(1, 3, 3, 3, 5), x = c(1, 2, 3, 4, 5), N = 5, xbar = 3)

The next step is to provide initial values for the MCMC. You should provide initial values for all random variables in the model. If you do not, BUGS will attempt to derive initial values for any random variables that are not initialised by using what you have provided. However sometimes it is unable to do so. The initial values are passed as a function with no arguments containing a named list with names corresponding to each random variable in the MCMC. R2OpenBUGS will call this function for every MCMC chain that is specified, generating initial values for each chain. Here we have integer values for alpha and beta. As we are running only 1 MCMC chain this is acceptable, but if we were running >1 chain, each chain would have the same starting values.

lineinits <- function() {    list(alpha = 1, beta = 1, tau = 1)}

Finally, let's run BUGS using the bugs(…) command and assign the output to an object called “lineout”. Note that we need to specify which nodes/parameters (random variables) in the model to save at each step of the MCMC. Basically, anything except data can be monitored by BUGS - model parameters, missing data values, functions of model parameters… This is a strength of BUGS since all of the parameters will have their posterior distribution samples saved for summarising, graphing etc. It can also place an overhead in data handling however: saving 1000 MCMC samples x k chains for elements such as fitted values e.g. mu[j] above can lead to a lot of data values. In this case it is trivial, but for other, larger datasets it can lead to problems.

Note also that here we only specify the number of MCMC iterations to be performed. By default R2OpenBUGS assumes that half of these iterations will be used as the burn-in. So in this instance the posterior distribution will be based on a sample of 500 MCMC iterations. In the case of this linear model that is probably sufficient - but also the number of iterations used as burn-in here is probably overkill.

## Uses default settings for n.burnin = n.iter/2; n.thin=1;
lineout <- bugs(data = linedata, inits = lineinits, parameters.to.save = c("alpha",     "beta", "sigma"), model.file = linemodel, n.chains = 1, n.iter = 10000)