MCMC.commonvariance = function(y,alpha,beta,lambda0,C0,theta0,V0,p,lambda,sigma2,burnin,M){ n = length(y) T = rep(0,n) mT = rep(0,n) chain = NULL for (iter in 1:(burnin+M)){ # sampling T's prob1 = p*dnorm(y,lambda[1],sqrt(sigma2)) prob2 = (1-p)*dnorm(y,lambda[2],sqrt(sigma2)) sprob = prob1+prob2 prob1 = prob1/sprob for (i in 1:n) T[i] = sample(1:2,size=1,prob=c(prob1[i],1-prob1[i])) # sampling sigma2 sigma2 = 1/rgamma(1,alpha+n/2,beta+sum((y-lambda[T])^2)/2) # sampling lambda[1] n1 = sum(T==1) n2 = n-n1 var = 1/(1/C0+n1/sigma2) mean = var*(lambda0/C0+sum(y[T==1])/sigma2) lambda[1] = rnorm(1,mean,sqrt(var)) # sampling p draw = rbeta(1,1+n1,1+n2) if (draw>0.5){ p = draw } #p = rbeta(1,1+n1,1+n2) # sampling theta var = 1/(1/V0+n2/sigma2) mean = var*(theta0/V0 + sum(y[T==2]-lambda[1])/sigma2) draw = rnorm(1,mean,sqrt(var)) if (draw>0){ theta = draw } lambda[2] = lambda[1]+theta if (iter>burnin){ chain = rbind(chain,c(p,lambda[1],theta,sigma2)) mT = mT + (T-1)/M } } return(list(p=chain[,1],lambda=chain[,2],theta=chain[,3],sigma2=chain[,4],mT=mT)) } MCMC.twovariances = function(y,alpha,beta,lambda0,C0,theta0,V0,p,lambda,sigma2,burnin,M){ n = length(y) T = rep(0,n) mT = rep(0,n) chain = NULL for (iter in 1:(burnin+M)){ # sampling T's prob1 = p*dnorm(y,lambda[1],sqrt(sigma2[1])) prob2 = (1-p)*dnorm(y,lambda[2],sqrt(sigma2[2])) sprob = prob1+prob2 prob1 = prob1/sprob for (i in 1:n) T[i] = sample(1:2,size=1,prob=c(prob1[i],1-prob1[i])) n1 = sum(T==1) n2 = n-n1 # sampling sigma2 sigma2[1] = 1/rgamma(1,alpha[1]+n1/2,beta[1]+sum((y[T==1]-lambda[1])^2)/2) sigma2[2] = 1/rgamma(1,alpha[2]+n2/2,beta[2]+sum((y[T==2]-lambda[2])^2)/2) # sampling lambda[1] var = 1/(1/C0+n1/sigma2[1]) mean = var*(lambda0/C0+sum(y[T==1])/sigma2[1]) lambda[1] = rnorm(1,mean,sqrt(var)) # sampling p draw = rbeta(1,1+n1,1+n2) if (draw>0.5){ p = draw } #p = rbeta(1,1+n1,1+n2) # sampling theta var = 1/(1/V0+n2/sigma2[2]) mean = var*(theta0/V0 + sum(y[T==2]-lambda[1])/sigma2[2]) draw = rnorm(1,mean,sqrt(var)) if (draw>0){ theta = draw } lambda[2] = lambda[1]+theta if (iter>burnin){ chain = rbind(chain,c(p,lambda[1],theta,sigma2[1],sigma2[2])) mT = mT + (T-1)/M } } return(list(p=chain[,1],lambda=chain[,2],theta=chain[,3],sigma21=chain[,4],sigma22=chain[,5],mT=mT)) } MCMC = function(y,alpha,beta,gamma,delta,p,lambda,k,sigma2,burnin,M){ T = rep(0,n) chain = NULL sigma2s = c(sigma2,k*sigma2) for (iter in 1:(burnin+M)){ # sampling T's prob1 = p*dnorm(y,lambda,sqrt(sigma2s[1])) prob2 = (1-p)*dnorm(y,lambda,sqrt(sigma2s[2])) sprob = prob1+prob2 prob1 = prob1/sprob for (i in 1:n) T[i] = sample(1:2,size=1,prob=c(prob1[i],1-prob1[i])) n1 = sum(T==1) ks = rep(1,n) ks[T==2] = k # sampling p p1 = rbeta(1,1+n1,1+n-n1) if (p1>0.5){ p = p1 } # sampling lambda var = 1/sum(1/sigma2s[T]) lambda = rnorm(1,var*sum(y/sigma2s[T]),sqrt(var)) # sampling sigma2 sigma2 = 1/rgamma(1,alpha+n/2,beta+sum((y-lambda)^2/ks)/2) sigma2s = c(sigma2,k*sigma2) # sampling k k1 = 1/rgamma(1,(n-n1-1)/2,sum((y[T==2]-lambda)^2)/(2*sigma2)) if (k1>1){ accept = min(1,(gamma-1)*log(k1-1)-delta*(k1-1)-(gamma-1)*log(k-1)+delta*(k-1)) if (log(runif(1))burnin){ chain = rbind(chain,c(p,lambda,sigma2,k)) } } return(chain) }