if (FALSE) {
## data generating
set.seed(1977)
Q <- 3
true.beta1 <- c(1, 1, 1) ; true.beta2 <- c(-1, -1, -1)
true.sigma2 <- c(2, 5); true.D1 <- diag(.5, Q); true.D2 <- diag(2.5, Q)
N=30; T=100;
NT <- N*T
x1 <- runif(NT, 1, 2)
x2 <- runif(NT, 1, 2)
X <- cbind(1, x1, x2); W <- X; y <- rep(NA, NT)
## true break numbers are one and at the center
break.point = rep(T/2, N); break.sigma=c(rep(1, N));
break.list <- rep(1, N)
id <- rep(1:N, each=NT/N)
K <- ncol(X);
ruler <- c(1:T)
## compute the weight for the break
W.mat <- matrix(NA, T, N)
for (i in 1:N){
W.mat[, i] <- pnorm((ruler-break.point[i])/break.sigma[i])
}
Weight <- as.vector(W.mat)
## data generating by weighting two means and variances
j = 1
for (i in 1:N){
Xi <- X[j:(j+T-1), ]
Wi <- W[j:(j+T-1), ]
true.V1 <- true.sigma2[1]*diag(T) + Wi%*%true.D1%*%t(Wi)
true.V2 <- true.sigma2[2]*diag(T) + Wi%*%true.D2%*%t(Wi)
true.mean1 <- Xi%*%true.beta1
true.mean2 <- Xi%*%true.beta2
weight <- Weight[j:(j+T-1)]
y[j:(j+T-1)] <- (1-weight)*true.mean1 + (1-weight)*chol(true.V1)%*%rnorm(T) +
weight*true.mean2 + weight*chol(true.V2)%*%rnorm(T)
j <- j + T
}
## model fitting
subject.id <- c(rep(1:N, each=T))
time.id <- c(rep(1:T, N))
## model fitting
G <- 100
b0 <- rep(0, K) ; B0 <- solve(diag(100, K))
c0 <- 2; d0 <- 2
r0 <- 5; R0 <- diag(c(1, 0.1, 0.1))
subject.id <- c(rep(1:N, each=T))
time.id <- c(rep(1:T, N))
out1 <- HMMpanelRE(subject.id, time.id, y, X, W, m=1,
mcmc=G, burnin=G, thin=1, verbose=G,
b0=b0, B0=B0, c0=c0, d0=d0, r0=r0, R0=R0)
## latent state changes
plotState(out1)
## print mcmc output
summary(out1)
}
Run the code above in your browser using DataLab