Skip to content

Instantly share code, notes, and snippets.

@mattblackwell
Created September 8, 2011 16:14
Show Gist options
  • Save mattblackwell/1203787 to your computer and use it in GitHub Desktop.
Save mattblackwell/1203787 to your computer and use it in GitHub Desktop.
Covariates in the Covariance Gibbs
## Function to implment the MCMC sampler
polling.gibbs.kal <- function(y, size, dates, states, big, ## data
F, X=NULL, XL = NULL, Z=NULL, ## covariates
m.0, C.0, b = NULL, B = NULL, ## hyper-parameters
c.0, d.0, e.0, f.0, g.0, h.0,
starts, phi, sigma.a, nu, beta, ## starting values
tune, A, tune.phi,
iters = 10000, thin = 1, burnin = 0) { ## mcmc stuff
aa <- proc.time()
## pre-gibbs cleaning (calculate variances, setup times, etc)
require(MCMCpack)
require(mvtnorm)
days <- dates- min(dates) + 1
T <- max(days)
num.states <- length(unique(states))
state.names <- unique(states)
## starting values
## gammas is the matrix (time by state) of state effects
gammas <- matrix(starts, nrow = T, ncol = num.states)
## Create "projection"
S <- exp(-mh(Z, diag(phi, ncol(Z))))
Delta <- sigma.a * S + diag(rep(nu, ncol(S)))
## Calculate number of draws
draws <- (iters-burnin)/thin
## set up acceptance ratios
acc.a <- acc.nu <- 0
acc.phi <- rep(0, length(phi))
## holder matrices
gamma.hold <- array(NA, dim=c(T,num.states,draws))
phi.hold <- matrix(NA, nrow = ncol(Z), ncol = draws)
sigma.a.hold <- rep(NA, times = draws)
nu.hold <- rep(NA, times = draws)
if (!is.null(X)) {
beta.hold <- matrix(NA, nrow = ncol(XL), ncol = draws)
}
for (i in 1:iters) {
## update alpha.0, gamma.s0
a <- matrix(0, nrow = T, ncol = num.states)
m <- matrix(0, nrow = T, ncol = num.states)
C <- array(0, dim = c(num.states,num.states,T))
R <- array(0, dim = c(num.states,num.states,T))
m[1,] <- m.0
C[,,1] <- C.0
## Begin Forward Filtering
for (t in 2:T) {
## y[,t] = F[[t]] %*% gammas[,t] + Sigma.t
## F relates the states (gammas) to the observations
## so we only need the rows that correspond to the
## states actually being used today
a[t,] <- m[t-1,]
R[,,t] <- C[,,t-1] + Delta
## slight modification for missing data.
## if there are no observations, then the posterior
## mean/variance is the prior mean/variance
if (nrow(big[[t]]) > 0) {
f <- F[[t]] %*% a[t,]
if (!is.null(X)) {
e <- big[[t]]$y - X[[t]] %*% beta - f
} else {
e <- big[[t]]$y - f
}
Sigma.t <- diag(big[[t]]$sigma.sq, nrow(big[[t]]))
Qt <- F[[t]] %*% R[, ,t] %*% t(F[[t]]) + Sigma.t
At <- R[,,t] %*% t(F[[t]]) %*% solve(Qt)
m[t,] <- a[t,] + At %*% e
C[,,t] <- R[,,t] - R[,,t] %*% t(F[[t]]) %*% solve(Qt) %*% F[[t]] %*% R[,,t]
} else {
m[t,] <- a[t,]
C[,,t] <- R[,,t]
}
}
## Begin Backward Sampling
gammas[T,] <- rmvnorm(1, m[T,], C[,,T])
if (!is.null(X)) {
beta.resid <- big[[T]]$y - F[[T]] %*% gammas[T,]
}
for (t in (T-1):1) {
B.t <- C[,,t] %*% solve(R[,,t+1])
h.t <- m[t,] + B.t %*% (gammas[t+1,] - a[t+1,])
H.t <- C[,,t] - B.t %*% R[,,t+1] %*% t(B.t)
gammas[t,] <- rmvnorm(1, h.t, H.t)
if (!is.null(X)) {
beta.resid <- c(beta.resid, big[[t]]$y - F[[t]] %*% gammas[t,])
}
}
if (!is.null(X)) {
## update betas (polling house effects)
Sigma.beta.star <- solve(t(XL) %*% diag(1/sigma.sq) %*% XL + B)
beta.star <- Sigma.beta.star %*% (t(XL) %*% diag(1/sigma.sq) %*% beta.resid + B %*% b)
beta <- mvrnorm(1, beta.star, Sigma.beta.star)
## normalize the house effects
beta <- beta - mean(beta)
}
## MH step for sigma.a
sigma.a.can <- rlnorm(1, log(sigma.a), sd = tune[1,1])
Delta.can <- sigma.a.can * S + diag(rep(nu, ncol(S)))
hold.can <- -dlnorm(sigma.a.can, log(sigma.a), tune[1,1], log = TRUE)
hold.cur <- -dlnorm(sigma.a, log(sigma.a.can), tune[1,1], log = TRUE)
hold.can <- hold.can + log(dinvgamma(sigma.a.can, c.0, d.0))
hold.cur <- hold.cur + log(dinvgamma(sigma.a, c.0, d.0))
for (t in 2:T) {
hold.can <- hold.can + dmvnorm(gammas[t,], gammas[t-1,], Delta.can, log = TRUE)
hold.cur <- hold.cur + dmvnorm(gammas[t,], gammas[t-1,], Delta, log = TRUE)
}
r.a <- exp(hold.can - hold.cur)
if (runif(1) <= r.a) {
sigma.a <- sigma.a.can
acc.a <- acc.a + 1
}
Delta <- sigma.a * S + diag(rep(nu, ncol(S)))
## MH step for nu ##
nu.can <- rlnorm(1, log(nu), sd = tune[2,2])
Delta.can <- sigma.a * S + diag(rep(nu.can, ncol(S)))
hold.can <- -dlnorm(nu.can, log(nu), tune[2, 2], log = TRUE)
hold.cur <- -dlnorm(nu, log(nu.can), tune[2, 2], log = TRUE)
hold.can <- hold.can + log(dinvgamma(nu.can, g.0, h.0))
hold.cur <- hold.cur + log(dinvgamma(nu, g.0, h.0))
for (t in 2:T) {
hold.can <- hold.can + dmvnorm(gammas[t,], gammas[t-1,], Delta.can, log = TRUE)
hold.cur <- hold.cur + dmvnorm(gammas[t,], gammas[t-1,], Delta, log = TRUE)
}
r.a <- exp(hold.can - hold.cur)
if (runif(1) <= r.a) {
nu <- nu.can
acc.nu <- acc.nu + 1
}
Delta <- sigma.a * S + diag(rep(nu, ncol(S)))
## MH for phi
for (p in 1:length(phi)) {
newphi <- rlnorm(1, log(phi[p]), sd = tune[p + 2, p + 2])
phi.can <- phi
phi.can[p] <- newphi
S.can <- exp(-mh(Z, diag(c(phi.can))))
Delta.can <- sigma.a * S.can + diag(rep(nu, ncol(S)))
## Jumping disttribution discount ##
hold.can <- -dlnorm(newphi, log(phi[p]), tune[p + 2, p + 2], log = TRUE)
hold.cur <- -dlnorm(phi[p], log(newphi), tune[p + 2, p + 2], log = TRUE)
## Priors
hold.can <- hold.can + dhalfcauchy(newphi, A, log = TRUE)
hold.cur <- hold.cur + dhalfcauchy(phi[p], A, log = TRUE)
for (t in 2:T) {
hold.can <- hold.can + dmvnorm(gammas[t,], gammas[t-1,], Delta.can, log = TRUE)
hold.cur <- hold.cur + dmvnorm(gammas[t,], gammas[t-1,], Delta, log = TRUE)
}
## Update phi[p] and Delta
r.a <- exp(hold.can - hold.cur)
if (runif(1) <= r.a) {
phi <- phi.can
acc.phi[p] <- acc.phi[p] + 1
S <- S.can
}
Delta <- sigma.a * S + diag(rep(nu, ncol(S)))
}
acc <- c(acc.a, acc.nu, acc.phi)/i
if ( (((i-burnin) %% thin)==0) & (i >burnin)) {
gamma.hold[,,((i-burnin)/thin)] <- gammas
phi.hold[,((i-burnin)/thin)] <- phi
sigma.a.hold[((i-burnin)/thin)] <- sigma.a
nu.hold[((i-burnin)/thin)] <- nu
if (!is.null(X)) {
beta.hold[,((i-burnin)/thin)] <- beta
save(list = c("gamma.hold", "beta.hold", "phi.hold", "sigma.a.hold", "nu.hold", "acc"),
file = "polls-out.RData")
}
}
if ((i %% 10)==0) cat("\n")
if ((i %% 100)==0) cat("Minutes Run: ",((proc.time()-aa)/60)[3],"\n")
cat(i," ")
#cat("\nCandidate: ",sqrt(sigma.a.can), sqrt(nu.can), phi.can, "\n")
cat("New: ", sqrt(sigma.a), sqrt(nu), phi, "\n")
cat("Acceptance rate: ", acc, "\n\n")
flush.console()
}
if (!is.null(X)) {
out <- list(gammas=gamma.hold, betas = beta.hold, phi = phi.hold,
sigma.a = sigma.a.hold, nu = nu.hold,
acc = acc)
} else {
out <- list(gammas=gamma.hold, phi = phi.hold,
sigma.a = sigma.a.hold, nu = nu.hold,
acc = acc)
}
return(out)
}
## Mahalanobis distance matrix
mh <- function(x, cov) {
S <- matrix(NA, nrow(x), nrow(x))
for (i in 1:nrow(S)) {
for (j in 1:i) {
S[i,j] <- sqrt(mahalanobis(x[i,], x[j,], cov=cov))
}
}
S[upper.tri(S)] <- t(S)[upper.tri(S)]
return(S)
}
dhalfcauchy <- function(x, A, log = FALSE) {
out <- -A - log(1 + x/(A^2))
if (!log) out <- exp(out)
return(out)
}
## Fake Covariate Data
N <- 9
Z <- as.matrix(expand.grid(1:3, 1:3))/10
s.a <- 0.01^2
nu <- 0.005^2
phi <- 0.1
D <- s.a * exp(-mh(Z, diag(phi, ncol(Z))))
## Create data with gammas as the latent variable, y as the observations
## I assume the observational variance is known
T <- 10
F.t <- list()
big <- list()
gs <- matrix(NA, ncol = T, nrow = N)
y <- matrix(NA, ncol = T, nrow = N)
gs[,1] <- runif(N, 0.4, 0.6)
y[,1] <- gs[,1] + rnorm(N, 0, sqrt((gs[,1] *(1-gs[,1]))/500))
for (t in 2:T) {
gs[,t] <- rmvnorm(1, gs[,t-1], D)
y[,t] <- gs[,t] + rnorm(N, 0, sqrt((gs[,t] *(1-gs[,t]))/500))
F.t[[t]] <- diag(N)
big[[t]] <- data.frame(y=y[,t], state=1:N, sigma.sq= sqrt(y[,t] * (1 - y[,t]))/sqrt(500))
}
y <- c(y)
dts <- rep(1:T, each = N)
sts <- rep(1:N, times = T)
## first day priors
m.0 <- c(rep(.5, times = N))
C.0 <- diag(c(rep(.1^2, times = N)))
## house effects
b <- rep(0, times = ncol(pollster.contrast))
B <- diag(rep(1/(.1^2), times = ncol(pollster.contrast)))
## common covariance
sigma.a <- 0.01^2
c.0 <- 1
d.0 <- 0.0005
tune.a <- 0.2
phi <- rep(.1, times = ncol(Z))
e.0 <- 1
f.0 <- 0.01
tune.phi <- c(0.5, 0.8, 0.8)##rep(0.5, ncol(Z))
## nugget
nu <- 0.005^2
g.0 <- 1
h.0 <- 0.0005
tune.nu <- .25
tune <- diag(c(tune.a, tune.nu, tune.phi))
out <- polling.gibbs.kal(y = y, dates = dts, big = big,states = sts, F = F.t, Z = Z,
m.0 =m.0, c.0 = c.0, d.0 = d.0, e.0 = e.0, f.0 = f.0,
C.0 = C.0, g.0 = g.0, h.0 = h.0,
starts = 0, phi = phi, sigma.a = sigma.a,
beta = NULL, nu = nu,
tune = tune, tune.phi = tune.phi, A = 5,
iters = 50, burnin=0,thin=1)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment