Commit df7868c3 authored by Chris Jewell's avatar Chris Jewell
Browse files

Merge branch 'modelselection' of fhm-chicas-code.lancs.ac.uk:jewell/zipBayes into modelselection

parents 427add33 56f7c125
......@@ -4,6 +4,25 @@
# of this file for data processing and running of
# MCMC.
ZipLik <- function(y, X, beta, gamma) {
# Split data
X.z <- X[y == 0, ]
Y.nz <- y[y != 0]
X.nz <- X[y != 0, ]
function() {
pi.z <- exp(X.z %*% gamma$val)
pi.z <- pi.z / (1 + pi.z)
li.z <- exp(X.z %*% beta$val)
pi.nz <- exp(X.nz %*% gamma$val)
pi.nz <- pi.nz / (1 + pi.nz)
li.nz <- exp(X.nz %*% beta$val)
zeropart <- sum(log(pi.z + (1 - pi.z) * exp(-li.z)))
nonzeropart <- sum(Y.nz * log(li.nz) - li.nz + log(1 - pi.nz))
zeropart + nonzeropart
}
}
####################
## MCMC algorithm ##
####################
......@@ -45,31 +64,33 @@ zipMCMC <-
# Parameters
beta <- new.env()
gamma <- new.env()
# Variable indicator -- start with full model.
z.beta <- rep(TRUE, ncol(X))
z.gamma <- rep(TRUE, ncol(X))
select.beta <- select.gamma <- rep(TRUE, ncol(X))
select.beta[force.beta] <- FALSE
select.gamma[force.gamma] <- FALSE
# Starting values
if(is.null(start.beta)) {
beta$val <- rep(0, ncol(X))
z.beta[select.beta] <- FALSE # Begin with null model
beta$val[1] <- glm(Y ~ 1, family = "poisson")$coef[1]
} else {
if(length(start.beta != ncol(X))) stop("Number of elements in start.beta not equal to columns in X")
if(length(start.beta) != ncol(X)) stop("Number of elements in start.beta not equal to columns in X")
beta$val <- start.beta
z.beta[beta$val==0] <- FALSE
}
if(is.null(start.gamma)) {
gamma$val <- rep(0, ncol(X))
z.gamma[select.gamma] <- FALSE # Begin with null model
gamma$val[1] <- glm((Y == 0) ~ 1, family = "binomial")$coef[1]
} else {
if(length(start.gamma != ncol(X))) stop("Number of elements in start.gamma not equal to columns in X")
if(length(start.gamma) != ncol(X)) stop("Number of elements in start.gamma not equal to columns in X")
gamma$val <- start.gamma
z.gamma[gamma$val==0] <- FALSE
}
# Priors
......@@ -80,24 +101,7 @@ zipMCMC <-
function()
sum(dnorm(val[val!=0], prior.gamma[val!=0, 1], prior.gamma[val!=0, 2], log = T)))
# Split data
X.z <- X[Y == 0, ]
Y.nz <- Y[Y != 0]
X.nz <- X[Y != 0, ]
llik <- function() {
g <- gamma$val
b <- beta$val
pi.z <- exp(X.z %*% g)
pi.z <- pi.z / (1 + pi.z)
li.z <- exp(X.z %*% b)
pi.nz <- exp(X.nz %*% g)
pi.nz <- pi.nz / (1 + pi.nz)
li.nz <- exp(X.nz %*% b)
zeropart <- sum(log(pi.z + (1 - pi.z) * exp(-li.z)))
nonzeropart <- sum(Y.nz * log(li.nz) - li.nz + log(1 - pi.nz))
zeropart + nonzeropart
}
llik <- ZipLik(Y, X, beta, gamma)
# Adaptive tuning
beta$tune <- rep(tune.beta, length(beta$val))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment