Nothing
R/supportings.R
In bcpmeta: Bayesian Multiple Changepoint Detection Using Metadata
Defines functions
loc2dirac
matA
lMarLik.normAR
lMarLik.EB.normAR
up.eta.MH.normAR
up.eta.RS.normAR
up.parameters.EB.normAR
lMarLik.phi.normAR
up.parameters.normAR
###################################################################
###### Supporting functions #######
###################################################################
###########################################################
## 1. Baic functions: marginal likelihood calculations ##
###########################################################
#### create meta or eta from non-zero-one sequences
loc2dirac = function(loc, n){
dirac = rep(0, n);
dirac[loc] = 1; dirac[1] = 0;
return(dirac);
}
#### create a (m + 1) dimensional squared matrix A, m != 0
matA = function(phi, m){
tmp = matrix(0, nrow = m + 1, ncol = m + 1);
indexoff = cbind(c(1:m, 2:(m + 1)), c(2:(m + 1), 1:m));
tmp[indexoff] = -phi;
diag(tmp) = c(1, rep(1 + phi^2, m - 1), 1);
return(tmp);
}
# Default parameter choices
# nu0 = 5; phi.upper = 0.99; phi.lower = -0.99;
# trend = TRUE or FALSE
lMarLik.normAR = function(Xtilde, eta, nu0, phi.lower, phi.upper, trend){
n = length(Xtilde);
m = sum(eta);
tau = which(eta == 1);
#### integrand p(Xtilde | phi, eta) * p(phi) ####
## these are fixed: eta, Xtilde, nu0
joint.Xtilde.phi = function(phi){
## calculate DtAD
if(m == 0)
DtAD = matrix((1 - phi^2) + (1 - phi)^2 * (n - 1), nrow = 1, ncol = 1);
if(m > 0){
DtAD = matA(phi, m);
diag(DtAD) = diag(DtAD) + (c(tau, n + 1) - c(1, tau) - 1) * (1- phi)^2;
}
## calculate XtdtAXtd and XtdtAT
predXtd = c(sqrt(1 - phi^2) * Xtilde[1], Xtilde[-1] - Xtilde[-n] * phi);
XtdtAXtd = sum(predXtd^2);
if(trend == TRUE){
predT = c(sqrt(1 - phi^2) * 1, (2:n) - (1:(n-1)) * phi);
TtAT = sum(predT^2);
XtdtAT = sum(predXtd * predT);
}
## calculate XtdtAD and TtAD
if(m == 0){
XtdtAD = sum(predXtd * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
XtdtAD = rep(0, m + 1);
XtdtAD[1] = sum(predXtd[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum1 = function(tauj){
return(
sum(predXtd[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
XtdtAD[2:m] = apply(tmp.index, 1, sum1);
}
XtdtAD[m + 1] = sum(predXtd[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
if(trend == TRUE){
if(m == 0){
TtAD = sum(predT * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
TtAD = rep(0, m + 1);
TtAD[1] = sum(predT[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum2 = function(tauj){
return(
sum(predT[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
TtAD[2:m] = apply(tmp.index, 1, sum2);
}
TtAD[m + 1] = sum(predT[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
}
invDtADnu0 = solve(DtAD + diag(1 / nu0, m + 1));
XtdtBXtd = XtdtAXtd - t(XtdtAD) %*% invDtADnu0 %*% XtdtAD;
if(trend == TRUE){
TtBT = TtAT - t(TtAD) %*% invDtADnu0 %*% TtAD;
XtdtBT = XtdtAT - t(XtdtAD) %*% invDtADnu0 %*% TtAD;
ljoint = 0.5 * (log(1 - phi^2) - log(det(DtAD + diag(1 / nu0, m + 1))) - log(TtBT)) - (n - 1) /2 * log(2 * pi) - (m + 1) / 2 * log(nu0) + lgamma((n - 1) / 2) - (n - 1) / 2 * log( (XtdtBXtd - XtdtBT^2 / TtBT) / 2 );
}
if(trend == FALSE){
ljoint = 0.5 * (log(1 - phi^2) - log(det(DtAD + diag(1 / nu0, m + 1)))) - n /2 * log(2 * pi) - (m + 1) / 2 * log(nu0) + lgamma(n / 2) - n / 2 * log( XtdtBXtd / 2 );
}
return(exp(ljoint));
}
lml = integrate(Vectorize(joint.Xtilde.phi), lower = phi.lower, upper = phi.upper);
return(list(lml = log(lml$value), lerr = log(lml$abs.error)));
}
########## Empirical Bayes (EB) estimates of sigmasq and phi ############
# Default parameter choices
# nu0 = 5; phi.upper = 0.99; phi.lower = -0.99;
# trend = TRUE or FALSE
lMarLik.EB.normAR = function(Xtilde, eta, nu0, phi.lower, phi.upper, trend){
n = length(Xtilde);
m = sum(eta);
tau = which(eta == 1);
#### integrand p(S | phi, gamma) * p(phi) ####
## these are fixed: gamma, S (which should be tilde(S) = S - mu0), nu0
joint.Xtilde.phi = function(phi){
## calculate DtAD
if(m == 0)
DtAD = matrix((1 - phi^2) + (1 - phi)^2 * (n - 1), nrow = 1, ncol = 1);
if(m > 0){
DtAD = matA(phi, m);
diag(DtAD) = diag(DtAD) + (c(tau, n + 1) - c(1, tau) - 1) * (1- phi)^2;
}
## calculate XtdtAXtd and XtdtAT
predXtd = c(sqrt(1 - phi^2) * Xtilde[1], Xtilde[-1] - Xtilde[-n] * phi);
XtdtAXtd = sum(predXtd^2);
if(trend == TRUE){
predT = c(sqrt(1 - phi^2) * 1, (2:n) - (1:(n-1)) * phi);
TtAT = sum(predT^2);
XtdtAT = sum(predXtd * predT);
}
## calculate XtdtAD and TtAD
if(m == 0){
XtdtAD = sum(predXtd * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
XtdtAD = rep(0, m + 1);
XtdtAD[1] = sum(predXtd[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum1 = function(tauj){
return(
sum(predXtd[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
XtdtAD[2:m] = apply(tmp.index, 1, sum1);
}
XtdtAD[m + 1] = sum(predXtd[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
if(trend == TRUE){
if(m == 0){
TtAD = sum(predT * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
TtAD = rep(0, m + 1);
TtAD[1] = sum(predT[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum2 = function(tauj){
return(
sum(predT[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
TtAD[2:m] = apply(tmp.index, 1, sum2);
}
TtAD[m + 1] = sum(predT[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
}
invDtADnu0 = solve(DtAD + diag(1 / nu0, m + 1));
XtdtBXtd = XtdtAXtd - t(XtdtAD) %*% invDtADnu0 %*% XtdtAD;
if(trend == TRUE){
TtBT = TtAT - t(TtAD) %*% invDtADnu0 %*% TtAD;
XtdtBT = XtdtAT - t(XtdtAD) %*% invDtADnu0 %*% TtAD;
ljoint = 0.5 * (log(1 - phi^2) - log(det(DtAD + diag(1 / nu0, m + 1))) - log(TtBT)) - (n - 1) /2 * log(2 * pi / (n - 1) * (XtdtBXtd - XtdtBT^2 / TtBT)) - (m + 1) / 2 * log(nu0) - (n - 1) / 2;
}
if(trend == FALSE){
ljoint = 0.5 * (log(1 - phi^2) - log(det(DtAD + diag(1 / nu0, m + 1)))) - n /2 * log(2 * pi / n * XtdtBXtd) - (m + 1) / 2 * log(nu0) - n / 2;
}
return(ljoint);
}
eb = optimize(joint.Xtilde.phi, interval = c(phi.lower, phi.upper), maximum = TRUE);
return(list(lml = c(eb$objective), phi.eb = c(eb$maximum)));
}
###################################################################
## 2. MCMC single iteration: Metropolis-Hastings update of eta ##
###################################################################
########## Update one random component of eta at a time ##########
## current = list(eta, lml, lpost, lerr, change.eta);
## if EB == TRUE, lerr is replaced by phi = phi.eb
up.eta.MH.normAR = function(Xtilde, meta, current, a1, a2, b1, b2, nu0, phi.lower, phi.upper, trend, EB){
eta = current$eta;
## randomly select the component to be flipped
n = length(Xtilde);
j = sample(2:n, 1);
# j = 51
eta2 = eta;
eta2[j] = 1 - eta[j];
m2.out = sum(eta2[(which(meta == 0))] == 1);
m2.in = sum(eta2[(which(meta == 1))] == 1);
if(EB == FALSE)
lmarlik2 = lMarLik.normAR(Xtilde, eta2, nu0, phi.lower, phi.upper, trend);
if(EB == TRUE)
lmarlik2 = lMarLik.EB.normAR(Xtilde, eta2, nu0, phi.lower, phi.upper, trend);
lpost2 = lmarlik2$lml + lbeta(a1 + m2.out, b1 + n - 1 - sum(meta) - m2.out) + lbeta(a2 + m2.in, b2 + sum(meta) - m2.in) - lbeta(a1, b1) - lbeta(a2, b2);
loga = lpost2 - current$lpost;
logu = log(runif(1));
if(loga > logu){
current$eta = eta2;
current$lml = lmarlik2$lml;
current$lpost = lpost2;
current$change.eta = TRUE;
if(EB == FALSE)
current$lerr = lmarlik2$lerr;
if(EB == TRUE)
current$phi = lmarlik2$phi.eb;
}
if(loga <= logu)
current$change.eta = FALSE;
return(current);
}
############ Simple random swapping update of eta ###########
## current = list(eta, lml, lpost, lerr, change.eta);
## if EB == TRUE, lerr is replaced by phi = phi.eb
up.eta.RS.normAR = function(Xtilde, meta, current, a1, a2, b1, b2, nu0, phi.lower, phi.upper, trend, EB){
eta = current$eta;
n = length(Xtilde);
if(sum(eta) > 0 && sum(eta) < n - 1){
if(sum(eta) == 1)
i = which(eta == 1);
if(sum(eta) > 1)
i = sample(which(eta == 1), 1);
if(sum(eta) == n - 2)
j = which(eta == 0)[2];
if(sum(eta) < n - 2)
j = sample(which(eta == 0)[-1], 1);
## randomly select the component to be flipped
eta2 = eta;
eta2[i] = 1 - eta[i];
eta2[j] = 1 - eta[j];
m2.out = sum(eta2[(which(meta == 0))] == 1);
m2.in = sum(eta2[(which(meta == 1))] == 1);
if(EB == FALSE)
lmarlik2 = lMarLik.normAR(Xtilde, eta2, nu0, phi.lower, phi.upper, trend);
if(EB == TRUE)
lmarlik2 = lMarLik.EB.normAR(Xtilde, eta2, nu0, phi.lower, phi.upper, trend);
lpost2 = lmarlik2$lml + lbeta(a1 + m2.out, b1 + n - 1 - sum(meta) - m2.out) + lbeta(a2 + m2.in, b2 + sum(meta) - m2.in) - lbeta(a1, b1) - lbeta(a2, b2);
loga = lpost2 - current$lpost;
logu = log(runif(1));
if(loga > logu){
current$eta = eta2;
current$lml = lmarlik2$lml;
current$lpost = lpost2;
current$change.eta = TRUE;
if(EB == FALSE)
current$lerr = lmarlik2$lerr;
if(EB == TRUE)
current$phi = lmarlik2$phi.eb;
}
if(loga <= logu)
current$change.eta = FALSE;
}
else{
current$change.eta = FALSE;
}
return(current);
}
###########################################################################
## 3. MCMC: model parameters Mu, Alpha, Sigmasq and Phi (if EB == FALSE) ##
###########################################################################
#############################
############ EB #############
#############################
## find EB estimate for sigamsq, and
## update alpha (if trend == TRUE), and mu (length N) using Gibbs sampler (no need to thin!)
## current = list(phi, sigmasq, alpha, mu);
up.parameters.EB.normAR = function(Xtilde, eta, phi.eb, iter = 1e3, nu0, trend){
n = length(Xtilde);
m = sum(eta);
tau = which(eta == 1);
phi = phi.eb;
if(m == 0)
DtAD = matrix((1 - phi^2) + (1 - phi)^2 * (n - 1), nrow = 1, ncol = 1);
if(m > 0){
DtAD = matA(phi, m);
diag(DtAD) = diag(DtAD) + (c(tau, n + 1) - c(1, tau) - 1) * (1- phi)^2;
}
## calculate XtdtAXtd and XtdtAT
predXtd = c(sqrt(1 - phi^2) * Xtilde[1], Xtilde[-1] - Xtilde[-n] * phi);
XtdtAXtd = sum(predXtd^2);
if(trend == TRUE){
predT = c(sqrt(1 - phi^2) * 1, (2:n) - (1:(n-1)) * phi);
TtAT = sum(predT^2);
XtdtAT = sum(predXtd * predT);
}
## calculate XtdtAD and TtAD
if(m == 0){
XtdtAD = sum(predXtd * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
XtdtAD = rep(0, m + 1);
XtdtAD[1] = sum(predXtd[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum1 = function(tauj){
return(
sum(predXtd[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
XtdtAD[2:m] = apply(tmp.index, 1, sum1);
}
XtdtAD[m + 1] = sum(predXtd[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
if(trend == TRUE){
if(m == 0){
TtAD = sum(predT * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
TtAD = rep(0, m + 1);
TtAD[1] = sum(predT[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum2 = function(tauj){
return(
sum(predT[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
TtAD[2:m] = apply(tmp.index, 1, sum2);
}
TtAD[m + 1] = sum(predT[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
}
invDtADnu0 = solve(DtAD + diag(1 / nu0, m + 1));
XtdtBXtd = XtdtAXtd - t(XtdtAD) %*% invDtADnu0 %*% XtdtAD;
if(trend == TRUE){
TtBT = TtAT - t(TtAD) %*% invDtADnu0 %*% TtAD;
XtdtBT = XtdtAT - t(XtdtAD) %*% invDtADnu0 %*% TtAD;
sigmasq = c((XtdtBXtd - XtdtBT^2 / TtBT) / (n - 1));
Alpha = rnorm(iter, XtdtBT / TtBT, sqrt(sigmasq / TtBT));
Mu = matrix(NA, nrow = iter, ncol = m + 1);
for(i in 1:iter)
Mu[i, ] = rmvnorm(1, invDtADnu0 %*% (XtdtAD - Alpha[i] * TtAD), sigmasq * invDtADnu0);
#Mu = Mu[, cumsum(eta) + 1];
#colnames(Mu) = paste('T', 1:n, sep = '');
return(list(sigmasq = sigmasq, Alpha = Alpha, Mu = Mu));
}
if(trend == FALSE){
sigmasq = c(XtdtBXtd / n);
Mu = rmvnorm(iter, invDtADnu0 %*% XtdtAD, sigmasq * invDtADnu0);
#Mu = Mu[, cumsum(eta) + 1];
#colnames(Mu) = paste('T', 1:n, sep = '');
return(list(sigmasq = sigmasq, Mu = Mu));
}
}
#############################
############ FB #############
#############################
## Calculate conditional marginal likelihood given phi and eta
lMarLik.phi.normAR = function(Xtilde, phi, eta, nu0, trend){
n = length(Xtilde);
m = sum(eta);
tau = which(eta == 1);
## calculate XtAX
if(m == 0)
DtAD = matrix((1 - phi^2) + (1 - phi)^2 * (n - 1), nrow = 1, ncol = 1);
if(m > 0){
DtAD = matA(phi, m);
diag(DtAD) = diag(DtAD) + (c(tau, n + 1) - c(1, tau) - 1) * (1- phi)^2;
}
## calculate XtdtAXtd and XtdtAT
predXtd = c(sqrt(1 - phi^2) * Xtilde[1], Xtilde[-1] - Xtilde[-n] * phi);
XtdtAXtd = sum(predXtd^2);
if(trend == TRUE){
predT = c(sqrt(1 - phi^2) * 1, (2:n) - (1:(n-1)) * phi);
TtAT = sum(predT^2);
XtdtAT = sum(predXtd * predT);
}
## calculate XtdtAD and TtAD
if(m == 0){
XtdtAD = sum(predXtd * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
XtdtAD = rep(0, m + 1);
XtdtAD[1] = sum(predXtd[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum1 = function(tauj){
return(
sum(predXtd[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
XtdtAD[2:m] = apply(tmp.index, 1, sum1);
}
XtdtAD[m + 1] = sum(predXtd[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
if(trend == TRUE){
if(m == 0){
TtAD = sum(predT * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
TtAD = rep(0, m + 1);
TtAD[1] = sum(predT[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum2 = function(tauj){
return(
sum(predT[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
TtAD[2:m] = apply(tmp.index, 1, sum2);
}
TtAD[m + 1] = sum(predT[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
}
invDtADnu0 = solve(DtAD + diag(1 / nu0, m + 1));
XtdtBXtd = XtdtAXtd - t(XtdtAD) %*% invDtADnu0 %*% XtdtAD;
if(trend == TRUE){
TtBT = TtAT - t(TtAD) %*% invDtADnu0 %*% TtAD;
XtdtBT = XtdtAT - t(XtdtAD) %*% invDtADnu0 %*% TtAD;
lml.phi = 0.5 * (log(1 - phi^2) - log(det(DtAD + diag(1 / nu0, m + 1))) - log(TtBT)) - (n - 1) /2 * log(2 * pi) - (m + 1) / 2 * log(nu0) + lgamma((n - 1) / 2) - (n - 1) / 2 * log( (XtdtBXtd - XtdtBT^2 / TtBT) / 2 );
return(list(lml.phi = lml.phi, invDtADnu0 = invDtADnu0, XtdtAD = XtdtAD, TtAD = TtAD, XtdtBT = XtdtBT, TtBT = TtBT, XtdtBXtd = XtdtBXtd));
}
if(trend == FALSE){
lml.phi = 0.5 * (log(1 - phi^2) - log(det(DtAD + diag(1 / nu0, m + 1)))) - n /2 * log(2 * pi) - (m + 1) / 2 * log(nu0) + lgamma(n / 2) - n / 2 * log( XtdtBXtd / 2 );
return(list(lml.phi = lml.phi, invDtADnu0 = invDtADnu0, XtdtAD = XtdtAD, XtdtBXtd = XtdtBXtd));
}
}
##############################################
## update phi using Metropolis-Hastings algorithm
## update sigmasq, alpha (if trend == TRUE), and mu (length N) using Gibbs sampler
## current = list(lml.phi, phi, sigmasq, alpha, mu, change.phi);
up.parameters.normAR = function(Xtilde, meta, eta, current, nu0, phi.lower, phi.upper, trend, sd.xi){
n = length(Xtilde);
m = sum(eta);
tau = which(eta == 1);
phi = current$phi;
lml = current$lml.phi;
##### M-H update of phi
xi = log((phi - phi.lower) / (phi.upper - phi));
xi2 = rnorm(1, phi, sd.xi);
phi2 = (exp(xi2) * phi.upper + phi.lower) / (exp(xi2) + 1);
lml2 = lMarLik.phi.normAR(Xtilde, phi2, eta, nu0, trend);
loga = lml2$lml.phi - lml$lml.phi + log(phi2 - phi.lower) + log(phi.upper - phi2) - log(phi - phi.lower) - log(phi.upper - phi);
logu = log(runif(1));
if(loga > logu){
current$phi = phi = phi2;
current$lml.phi = lml = lml2;
current$change.phi = current$change.phi + 1;
}
#### Gibbs update of sigma, alpha, mu (if trend == TRUE)
if(trend == TRUE){
current$sigmasq = 1 / rgamma(1, 0.5 * (n - 1), 0.5 * (lml$XtdtBXtd - lml$XtdtBT^2 / lml$TtBT) );
current$alpha = rnorm(1, lml$XtdtBT / lml$TtBT, sqrt(current$sigmasq / lml$TtBT));
mu = rmvnorm(1, lml$invDtADnu0 %*% (lml$XtdtAD - current$alpha * lml$TtAD), current$sigmasq * lml$invDtADnu0);
current$mu = mu;
#current$mu = mu[cumsum(eta) + 1];
}
#### Gibbs update of sigma, alpha, mu (if trend == FALSE)
if(trend == FALSE){
current$sigmasq = 1 / rgamma(1, 0.5 * n, 0.5 * lml$XtdtBXtd);
mu = rmvnorm(1, lml$invDtADnu0 %*% lml$XtdtAD, current$sigmasq * lml$invDtADnu0);
current$mu = mu;
#current$mu = mu[cumsum(eta) + 1];
}
return( current );
}
Try the bcpmeta package in your browser
Any scripts or data that you put into this service are public.
bcpmeta documentation built on May 1, 2019, 6:48 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions loc2dirac matA lMarLik.normAR lMarLik.EB.normAR up.eta.MH.normAR up.eta.RS.normAR up.parameters.EB.normAR lMarLik.phi.normAR up.parameters.normAR
###################################################################
###### Supporting functions #######
###################################################################
###########################################################
## 1. Baic functions: marginal likelihood calculations ##
###########################################################
#### create meta or eta from non-zero-one sequences
loc2dirac = function(loc, n){
dirac = rep(0, n);
dirac[loc] = 1; dirac[1] = 0;
return(dirac);
}
#### create a (m + 1) dimensional squared matrix A, m != 0
matA = function(phi, m){
tmp = matrix(0, nrow = m + 1, ncol = m + 1);
indexoff = cbind(c(1:m, 2:(m + 1)), c(2:(m + 1), 1:m));
tmp[indexoff] = -phi;
diag(tmp) = c(1, rep(1 + phi^2, m - 1), 1);
return(tmp);
}
# Default parameter choices
# nu0 = 5; phi.upper = 0.99; phi.lower = -0.99;
# trend = TRUE or FALSE
lMarLik.normAR = function(Xtilde, eta, nu0, phi.lower, phi.upper, trend){
n = length(Xtilde);
m = sum(eta);
tau = which(eta == 1);
#### integrand p(Xtilde | phi, eta) * p(phi) ####
## these are fixed: eta, Xtilde, nu0
joint.Xtilde.phi = function(phi){
## calculate DtAD
if(m == 0)
DtAD = matrix((1 - phi^2) + (1 - phi)^2 * (n - 1), nrow = 1, ncol = 1);
if(m > 0){
DtAD = matA(phi, m);
diag(DtAD) = diag(DtAD) + (c(tau, n + 1) - c(1, tau) - 1) * (1- phi)^2;
}
## calculate XtdtAXtd and XtdtAT
predXtd = c(sqrt(1 - phi^2) * Xtilde[1], Xtilde[-1] - Xtilde[-n] * phi);
XtdtAXtd = sum(predXtd^2);
if(trend == TRUE){
predT = c(sqrt(1 - phi^2) * 1, (2:n) - (1:(n-1)) * phi);
TtAT = sum(predT^2);
XtdtAT = sum(predXtd * predT);
}
## calculate XtdtAD and TtAD
if(m == 0){
XtdtAD = sum(predXtd * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
XtdtAD = rep(0, m + 1);
XtdtAD[1] = sum(predXtd[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum1 = function(tauj){
return(
sum(predXtd[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
XtdtAD[2:m] = apply(tmp.index, 1, sum1);
}
XtdtAD[m + 1] = sum(predXtd[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
if(trend == TRUE){
if(m == 0){
TtAD = sum(predT * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
TtAD = rep(0, m + 1);
TtAD[1] = sum(predT[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum2 = function(tauj){
return(
sum(predT[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
TtAD[2:m] = apply(tmp.index, 1, sum2);
}
TtAD[m + 1] = sum(predT[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
}
invDtADnu0 = solve(DtAD + diag(1 / nu0, m + 1));
XtdtBXtd = XtdtAXtd - t(XtdtAD) %*% invDtADnu0 %*% XtdtAD;
if(trend == TRUE){
TtBT = TtAT - t(TtAD) %*% invDtADnu0 %*% TtAD;
XtdtBT = XtdtAT - t(XtdtAD) %*% invDtADnu0 %*% TtAD;
ljoint = 0.5 * (log(1 - phi^2) - log(det(DtAD + diag(1 / nu0, m + 1))) - log(TtBT)) - (n - 1) /2 * log(2 * pi) - (m + 1) / 2 * log(nu0) + lgamma((n - 1) / 2) - (n - 1) / 2 * log( (XtdtBXtd - XtdtBT^2 / TtBT) / 2 );
}
if(trend == FALSE){
ljoint = 0.5 * (log(1 - phi^2) - log(det(DtAD + diag(1 / nu0, m + 1)))) - n /2 * log(2 * pi) - (m + 1) / 2 * log(nu0) + lgamma(n / 2) - n / 2 * log( XtdtBXtd / 2 );
}
return(exp(ljoint));
}
lml = integrate(Vectorize(joint.Xtilde.phi), lower = phi.lower, upper = phi.upper);
return(list(lml = log(lml$value), lerr = log(lml$abs.error)));
}
########## Empirical Bayes (EB) estimates of sigmasq and phi ############
# Default parameter choices
# nu0 = 5; phi.upper = 0.99; phi.lower = -0.99;
# trend = TRUE or FALSE
lMarLik.EB.normAR = function(Xtilde, eta, nu0, phi.lower, phi.upper, trend){
n = length(Xtilde);
m = sum(eta);
tau = which(eta == 1);
#### integrand p(S | phi, gamma) * p(phi) ####
## these are fixed: gamma, S (which should be tilde(S) = S - mu0), nu0
joint.Xtilde.phi = function(phi){
## calculate DtAD
if(m == 0)
DtAD = matrix((1 - phi^2) + (1 - phi)^2 * (n - 1), nrow = 1, ncol = 1);
if(m > 0){
DtAD = matA(phi, m);
diag(DtAD) = diag(DtAD) + (c(tau, n + 1) - c(1, tau) - 1) * (1- phi)^2;
}
## calculate XtdtAXtd and XtdtAT
predXtd = c(sqrt(1 - phi^2) * Xtilde[1], Xtilde[-1] - Xtilde[-n] * phi);
XtdtAXtd = sum(predXtd^2);
if(trend == TRUE){
predT = c(sqrt(1 - phi^2) * 1, (2:n) - (1:(n-1)) * phi);
TtAT = sum(predT^2);
XtdtAT = sum(predXtd * predT);
}
## calculate XtdtAD and TtAD
if(m == 0){
XtdtAD = sum(predXtd * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
XtdtAD = rep(0, m + 1);
XtdtAD[1] = sum(predXtd[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum1 = function(tauj){
return(
sum(predXtd[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
XtdtAD[2:m] = apply(tmp.index, 1, sum1);
}
XtdtAD[m + 1] = sum(predXtd[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
if(trend == TRUE){
if(m == 0){
TtAD = sum(predT * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
TtAD = rep(0, m + 1);
TtAD[1] = sum(predT[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum2 = function(tauj){
return(
sum(predT[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
TtAD[2:m] = apply(tmp.index, 1, sum2);
}
TtAD[m + 1] = sum(predT[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
}
invDtADnu0 = solve(DtAD + diag(1 / nu0, m + 1));
XtdtBXtd = XtdtAXtd - t(XtdtAD) %*% invDtADnu0 %*% XtdtAD;
if(trend == TRUE){
TtBT = TtAT - t(TtAD) %*% invDtADnu0 %*% TtAD;
XtdtBT = XtdtAT - t(XtdtAD) %*% invDtADnu0 %*% TtAD;
ljoint = 0.5 * (log(1 - phi^2) - log(det(DtAD + diag(1 / nu0, m + 1))) - log(TtBT)) - (n - 1) /2 * log(2 * pi / (n - 1) * (XtdtBXtd - XtdtBT^2 / TtBT)) - (m + 1) / 2 * log(nu0) - (n - 1) / 2;
}
if(trend == FALSE){
ljoint = 0.5 * (log(1 - phi^2) - log(det(DtAD + diag(1 / nu0, m + 1)))) - n /2 * log(2 * pi / n * XtdtBXtd) - (m + 1) / 2 * log(nu0) - n / 2;
}
return(ljoint);
}
eb = optimize(joint.Xtilde.phi, interval = c(phi.lower, phi.upper), maximum = TRUE);
return(list(lml = c(eb$objective), phi.eb = c(eb$maximum)));
}
###################################################################
## 2. MCMC single iteration: Metropolis-Hastings update of eta ##
###################################################################
########## Update one random component of eta at a time ##########
## current = list(eta, lml, lpost, lerr, change.eta);
## if EB == TRUE, lerr is replaced by phi = phi.eb
up.eta.MH.normAR = function(Xtilde, meta, current, a1, a2, b1, b2, nu0, phi.lower, phi.upper, trend, EB){
eta = current$eta;
## randomly select the component to be flipped
n = length(Xtilde);
j = sample(2:n, 1);
# j = 51
eta2 = eta;
eta2[j] = 1 - eta[j];
m2.out = sum(eta2[(which(meta == 0))] == 1);
m2.in = sum(eta2[(which(meta == 1))] == 1);
if(EB == FALSE)
lmarlik2 = lMarLik.normAR(Xtilde, eta2, nu0, phi.lower, phi.upper, trend);
if(EB == TRUE)
lmarlik2 = lMarLik.EB.normAR(Xtilde, eta2, nu0, phi.lower, phi.upper, trend);
lpost2 = lmarlik2$lml + lbeta(a1 + m2.out, b1 + n - 1 - sum(meta) - m2.out) + lbeta(a2 + m2.in, b2 + sum(meta) - m2.in) - lbeta(a1, b1) - lbeta(a2, b2);
loga = lpost2 - current$lpost;
logu = log(runif(1));
if(loga > logu){
current$eta = eta2;
current$lml = lmarlik2$lml;
current$lpost = lpost2;
current$change.eta = TRUE;
if(EB == FALSE)
current$lerr = lmarlik2$lerr;
if(EB == TRUE)
current$phi = lmarlik2$phi.eb;
}
if(loga <= logu)
current$change.eta = FALSE;
return(current);
}
############ Simple random swapping update of eta ###########
## current = list(eta, lml, lpost, lerr, change.eta);
## if EB == TRUE, lerr is replaced by phi = phi.eb
up.eta.RS.normAR = function(Xtilde, meta, current, a1, a2, b1, b2, nu0, phi.lower, phi.upper, trend, EB){
eta = current$eta;
n = length(Xtilde);
if(sum(eta) > 0 && sum(eta) < n - 1){
if(sum(eta) == 1)
i = which(eta == 1);
if(sum(eta) > 1)
i = sample(which(eta == 1), 1);
if(sum(eta) == n - 2)
j = which(eta == 0)[2];
if(sum(eta) < n - 2)
j = sample(which(eta == 0)[-1], 1);
## randomly select the component to be flipped
eta2 = eta;
eta2[i] = 1 - eta[i];
eta2[j] = 1 - eta[j];
m2.out = sum(eta2[(which(meta == 0))] == 1);
m2.in = sum(eta2[(which(meta == 1))] == 1);
if(EB == FALSE)
lmarlik2 = lMarLik.normAR(Xtilde, eta2, nu0, phi.lower, phi.upper, trend);
if(EB == TRUE)
lmarlik2 = lMarLik.EB.normAR(Xtilde, eta2, nu0, phi.lower, phi.upper, trend);
lpost2 = lmarlik2$lml + lbeta(a1 + m2.out, b1 + n - 1 - sum(meta) - m2.out) + lbeta(a2 + m2.in, b2 + sum(meta) - m2.in) - lbeta(a1, b1) - lbeta(a2, b2);
loga = lpost2 - current$lpost;
logu = log(runif(1));
if(loga > logu){
current$eta = eta2;
current$lml = lmarlik2$lml;
current$lpost = lpost2;
current$change.eta = TRUE;
if(EB == FALSE)
current$lerr = lmarlik2$lerr;
if(EB == TRUE)
current$phi = lmarlik2$phi.eb;
}
if(loga <= logu)
current$change.eta = FALSE;
}
else{
current$change.eta = FALSE;
}
return(current);
}
###########################################################################
## 3. MCMC: model parameters Mu, Alpha, Sigmasq and Phi (if EB == FALSE) ##
###########################################################################
#############################
############ EB #############
#############################
## find EB estimate for sigamsq, and
## update alpha (if trend == TRUE), and mu (length N) using Gibbs sampler (no need to thin!)
## current = list(phi, sigmasq, alpha, mu);
up.parameters.EB.normAR = function(Xtilde, eta, phi.eb, iter = 1e3, nu0, trend){
n = length(Xtilde);
m = sum(eta);
tau = which(eta == 1);
phi = phi.eb;
if(m == 0)
DtAD = matrix((1 - phi^2) + (1 - phi)^2 * (n - 1), nrow = 1, ncol = 1);
if(m > 0){
DtAD = matA(phi, m);
diag(DtAD) = diag(DtAD) + (c(tau, n + 1) - c(1, tau) - 1) * (1- phi)^2;
}
## calculate XtdtAXtd and XtdtAT
predXtd = c(sqrt(1 - phi^2) * Xtilde[1], Xtilde[-1] - Xtilde[-n] * phi);
XtdtAXtd = sum(predXtd^2);
if(trend == TRUE){
predT = c(sqrt(1 - phi^2) * 1, (2:n) - (1:(n-1)) * phi);
TtAT = sum(predT^2);
XtdtAT = sum(predXtd * predT);
}
## calculate XtdtAD and TtAD
if(m == 0){
XtdtAD = sum(predXtd * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
XtdtAD = rep(0, m + 1);
XtdtAD[1] = sum(predXtd[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum1 = function(tauj){
return(
sum(predXtd[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
XtdtAD[2:m] = apply(tmp.index, 1, sum1);
}
XtdtAD[m + 1] = sum(predXtd[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
if(trend == TRUE){
if(m == 0){
TtAD = sum(predT * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
TtAD = rep(0, m + 1);
TtAD[1] = sum(predT[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum2 = function(tauj){
return(
sum(predT[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
TtAD[2:m] = apply(tmp.index, 1, sum2);
}
TtAD[m + 1] = sum(predT[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
}
invDtADnu0 = solve(DtAD + diag(1 / nu0, m + 1));
XtdtBXtd = XtdtAXtd - t(XtdtAD) %*% invDtADnu0 %*% XtdtAD;
if(trend == TRUE){
TtBT = TtAT - t(TtAD) %*% invDtADnu0 %*% TtAD;
XtdtBT = XtdtAT - t(XtdtAD) %*% invDtADnu0 %*% TtAD;
sigmasq = c((XtdtBXtd - XtdtBT^2 / TtBT) / (n - 1));
Alpha = rnorm(iter, XtdtBT / TtBT, sqrt(sigmasq / TtBT));
Mu = matrix(NA, nrow = iter, ncol = m + 1);
for(i in 1:iter)
Mu[i, ] = rmvnorm(1, invDtADnu0 %*% (XtdtAD - Alpha[i] * TtAD), sigmasq * invDtADnu0);
#Mu = Mu[, cumsum(eta) + 1];
#colnames(Mu) = paste('T', 1:n, sep = '');
return(list(sigmasq = sigmasq, Alpha = Alpha, Mu = Mu));
}
if(trend == FALSE){
sigmasq = c(XtdtBXtd / n);
Mu = rmvnorm(iter, invDtADnu0 %*% XtdtAD, sigmasq * invDtADnu0);
#Mu = Mu[, cumsum(eta) + 1];
#colnames(Mu) = paste('T', 1:n, sep = '');
return(list(sigmasq = sigmasq, Mu = Mu));
}
}
#############################
############ FB #############
#############################
## Calculate conditional marginal likelihood given phi and eta
lMarLik.phi.normAR = function(Xtilde, phi, eta, nu0, trend){
n = length(Xtilde);
m = sum(eta);
tau = which(eta == 1);
## calculate XtAX
if(m == 0)
DtAD = matrix((1 - phi^2) + (1 - phi)^2 * (n - 1), nrow = 1, ncol = 1);
if(m > 0){
DtAD = matA(phi, m);
diag(DtAD) = diag(DtAD) + (c(tau, n + 1) - c(1, tau) - 1) * (1- phi)^2;
}
## calculate XtdtAXtd and XtdtAT
predXtd = c(sqrt(1 - phi^2) * Xtilde[1], Xtilde[-1] - Xtilde[-n] * phi);
XtdtAXtd = sum(predXtd^2);
if(trend == TRUE){
predT = c(sqrt(1 - phi^2) * 1, (2:n) - (1:(n-1)) * phi);
TtAT = sum(predT^2);
XtdtAT = sum(predXtd * predT);
}
## calculate XtdtAD and TtAD
if(m == 0){
XtdtAD = sum(predXtd * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
XtdtAD = rep(0, m + 1);
XtdtAD[1] = sum(predXtd[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum1 = function(tauj){
return(
sum(predXtd[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
XtdtAD[2:m] = apply(tmp.index, 1, sum1);
}
XtdtAD[m + 1] = sum(predXtd[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
if(trend == TRUE){
if(m == 0){
TtAD = sum(predT * c(sqrt(1 - phi^2), rep(1 - phi, n - 1)));
}
if(m > 0){
TtAD = rep(0, m + 1);
TtAD[1] = sum(predT[1 : tau[1]] * c(sqrt(1 - phi^2), rep(1 - phi, tau[1] - 2), -phi));
tmp.index = cbind(tau[-m], tau[-1]);
sum2 = function(tauj){
return(
sum(predT[tauj[1] : tauj[2]] * c(1, rep(1 - phi, tauj[2] - tauj[1] - 1), -phi))
);
}
if(length(tmp.index) > 0){
TtAD[2:m] = apply(tmp.index, 1, sum2);
}
TtAD[m + 1] = sum(predT[tau[m]:n] * c(1, rep(1 - phi, n - tau[m])));
}
}
invDtADnu0 = solve(DtAD + diag(1 / nu0, m + 1));
XtdtBXtd = XtdtAXtd - t(XtdtAD) %*% invDtADnu0 %*% XtdtAD;
if(trend == TRUE){
TtBT = TtAT - t(TtAD) %*% invDtADnu0 %*% TtAD;
XtdtBT = XtdtAT - t(XtdtAD) %*% invDtADnu0 %*% TtAD;
lml.phi = 0.5 * (log(1 - phi^2) - log(det(DtAD + diag(1 / nu0, m + 1))) - log(TtBT)) - (n - 1) /2 * log(2 * pi) - (m + 1) / 2 * log(nu0) + lgamma((n - 1) / 2) - (n - 1) / 2 * log( (XtdtBXtd - XtdtBT^2 / TtBT) / 2 );
return(list(lml.phi = lml.phi, invDtADnu0 = invDtADnu0, XtdtAD = XtdtAD, TtAD = TtAD, XtdtBT = XtdtBT, TtBT = TtBT, XtdtBXtd = XtdtBXtd));
}
if(trend == FALSE){
lml.phi = 0.5 * (log(1 - phi^2) - log(det(DtAD + diag(1 / nu0, m + 1)))) - n /2 * log(2 * pi) - (m + 1) / 2 * log(nu0) + lgamma(n / 2) - n / 2 * log( XtdtBXtd / 2 );
return(list(lml.phi = lml.phi, invDtADnu0 = invDtADnu0, XtdtAD = XtdtAD, XtdtBXtd = XtdtBXtd));
}
}
##############################################
## update phi using Metropolis-Hastings algorithm
## update sigmasq, alpha (if trend == TRUE), and mu (length N) using Gibbs sampler
## current = list(lml.phi, phi, sigmasq, alpha, mu, change.phi);
up.parameters.normAR = function(Xtilde, meta, eta, current, nu0, phi.lower, phi.upper, trend, sd.xi){
n = length(Xtilde);
m = sum(eta);
tau = which(eta == 1);
phi = current$phi;
lml = current$lml.phi;
##### M-H update of phi
xi = log((phi - phi.lower) / (phi.upper - phi));
xi2 = rnorm(1, phi, sd.xi);
phi2 = (exp(xi2) * phi.upper + phi.lower) / (exp(xi2) + 1);
lml2 = lMarLik.phi.normAR(Xtilde, phi2, eta, nu0, trend);
loga = lml2$lml.phi - lml$lml.phi + log(phi2 - phi.lower) + log(phi.upper - phi2) - log(phi - phi.lower) - log(phi.upper - phi);
logu = log(runif(1));
if(loga > logu){
current$phi = phi = phi2;
current$lml.phi = lml = lml2;
current$change.phi = current$change.phi + 1;
}
#### Gibbs update of sigma, alpha, mu (if trend == TRUE)
if(trend == TRUE){
current$sigmasq = 1 / rgamma(1, 0.5 * (n - 1), 0.5 * (lml$XtdtBXtd - lml$XtdtBT^2 / lml$TtBT) );
current$alpha = rnorm(1, lml$XtdtBT / lml$TtBT, sqrt(current$sigmasq / lml$TtBT));
mu = rmvnorm(1, lml$invDtADnu0 %*% (lml$XtdtAD - current$alpha * lml$TtAD), current$sigmasq * lml$invDtADnu0);
current$mu = mu;
#current$mu = mu[cumsum(eta) + 1];
}
#### Gibbs update of sigma, alpha, mu (if trend == FALSE)
if(trend == FALSE){
current$sigmasq = 1 / rgamma(1, 0.5 * n, 0.5 * lml$XtdtBXtd);
mu = rmvnorm(1, lml$invDtADnu0 %*% lml$XtdtAD, current$sigmasq * lml$invDtADnu0);
current$mu = mu;
#current$mu = mu[cumsum(eta) + 1];
}
return( current );
}
Try the bcpmeta package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.