I try to do a test for dirichlet process for Multivariate normal, but Winbugs
always says "expected multivariate node", does that mean I miss something at
initialization? I will really appreciate the help to solve this problem

Here is the R code, and Winbugs code.



model
{
        for(i in 1:N){
        y[i,1:2] ~ dmnorm(mu[i,],tau[i,,])      
        S[i] ~ dcat(pi[])
        mu[i,1:2] <- mu.star[S[i],]
        tau[i,1:2,1:2] <- tau.star[S[i],,]}

        # Constructive DPP
        
        # Stick breaking prior
        p[1] <- r[1]
        for (j in 2:C) {p[j] <- r[j]*(1-r[j-1])*p[j-1]/r[j-1]}
        p.sum <- sum(p[])               
        for (j in 1:C) {r[j] ~ dbeta(1,alpha); pi[j] <- p[j]/p.sum}
        
        # Baseline distribution
        for (j in 1:C) {mu.star[j,1:2] ~
dmnorm(theta[],tau.star[j,,]);tau.star[j,1:2,1:2] ~ dwish(T[,],3)}
        theta[1:2] ~ dmnorm(theta0[],S2[,])
        T[1:2,1:2] ~ dwish(S3[,],3)

        # DPP Precision Parameter               
        alpha ~ dgamma(1,1)

        # Programming for calculating summary statistics
        for(i in 1:N) {for (j in 1:C) {SC[i,j] <- equals(j,S[i])}}

        # total clusters K                      
        for (j in 1:C) {cl[j] <- step(sum(SC[,j])-1)}
        K <- sum(cl[])
}



library(R2WinBUGS)
library(MASS)
library(MCMCpack)
library(coda)

#Generate synthetic data
n=50
S=matrix(c(1,.2,.2,2),nrow=2)
y1=mvrnorm(n,c(0,0),S)
y2=mvrnorm(n,c(3,3),S)
y3=mvrnorm(n,c(-3,-3),S)
y=rbind(y1,y2,y3)

N=nrow(y)
C=6
theta0=as.vector(c(0,0))
S2=matrix(c(1,0,0,2),nrow=2)
S3=matrix(c(1,0,0,1),nrow=2)

data=list("y","N","C","theta0","S2","S3")
inits=function(){list(alpha=rgamma(1,1.5,1),theta=mvrnorm(1,theta0,S2),T=rwish(3,S3),K=6)}
dp_norm.sim=bugs(data,inits,model.file="dp_norm.txt",parameters=c("alpha","K"),
n.chains=1,n.iter=500,n.burnin=50,n.thin=1,bugs.directory="C:/Program
Files/WinBUGS14/",codaPkg=T,debug=T)





-- 
View this message in context: 
http://n4.nabble.com/a-small-question-about-R-with-Winbugs-tp1788668p1788668.html
Sent from the R help mailing list archive at Nabble.com.

______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to