Nothing
Log.post.theta <- function(K, od, theta) {
# theta: list, sigma, tau each is a K-vec, rho: K * K matrix
sigma2 <- theta[[1]] # vector
tau <- theta[[2]] # actually is inverse tau!!
rho <- theta[[3]] # K * K matrix: (1, rho, rho,1)
if(K > 1){
rho = matrix(rho, nrow = K, ncol =K)
diag(rho) = 1
}else rho = matrix(rho)
# sorting time given delta
time.comb <- od$time.comb
filter_ind <- od$filter_ind
var.lc.comb <- od$var
lc.comb <- od$lc # observations
leng.time = length(time.comb)
Sigma <- list() # T: K*K, sigma t|t-1
Sigma.star <- list() # T: K*K, sigma t|t
Sigma[[1]] = t(rho * sqrt(sigma2)) * sqrt(sigma2)
temp = matrix(0, K, K)
for(i in 1:K)
{
for(j in 1:K){
temp[i, j] = tau[i] + tau[j]
}
}
Sigma[[1]] = Sigma[[1]] /temp # Q
# x.hat.i, i = 1, 2, ..., 2n to be saved
mu.i <- matrix(0, leng.time, K) # T * K, mu t|t-1
mu.star.i <- matrix(0, leng.time, K) # mu t|t
# a.i, i = 2, ..., 2n
delta.t <- diff(time.comb)
a.i <- exp( -delta.t %*% t(tau)) # (T-1) * K
for (t in 1 : leng.time) {
if(t>1) Sigma[[t]] <- t(a.i[t - 1, ] *Sigma.star[[t-1]]) * a.i[t - 1, ] + Sigma[[1]] * (1 - exp(-temp*delta.t[t-1]))
k = ceiling(filter_ind[t]/2)
Sigma.star[[t]] <- Sigma[[t]] - Sigma[[t]][, k] %*% t(Sigma[[t]][,k]) /(Sigma[[t]][k, k] + var.lc.comb[t])
if(t>1) mu.i[t,] <- a.i[t - 1, ] * mu.star.i[t-1, ]
mu.star.i[t,] <- mu.i[t,] + Sigma[[t]][, k] /(Sigma[[t]][k, k] + var.lc.comb[t]) * (lc.comb[t] - mu.i[t,k])
}
# log-likelihood
mus <- sapply(1:leng.time, function(t) mu.i[t, ceiling(filter_ind[t]/2)])
vars <- sapply(1:leng.time, function(t) Sigma[[t]][ceiling(filter_ind[t]/2), ceiling(filter_ind[t]/2)] + var.lc.comb[t])
if(any(vars<0))
{
print(vars[vars<0])
print(time.comb[which(vars<0)])
}
return(sum(dnorm(lc.comb, mean = mus, sd = sqrt(vars), log = TRUE)))
}
postTheta.multi <- function (K, od, previous.theta,
tau.thresh, sigma2.thresh, rho.thresh,
tau.jump,sigma2.jump, rho.proposal.scale,
tau.prior.a, tau.prior.b, sigma.prior.a,
sigma.prior.b, rho.prior.a, rho.prior.b) {
# updating sigma
sigma2 = previous.theta[[1]]
sigma2.p <- exp(log(sigma2) + sigma2.jump)
for(k in 1:K)
{
theta = previous.theta
theta[[1]][k] = sigma2.p[k]
l.metrop <- Log.post.theta(K, od, theta) - Log.post.theta(K, od, previous.theta)
l.metrop <- l.metrop - sigma.prior.b/sigma2.p[k] - (sigma.prior.a + 1) * log(sigma2.p[k])
+ sigma.prior.b/sigma2[k] + (sigma.prior.a + 1) * log(sigma2[k])
l.hastings <- log(sigma2.p[k]) - log(sigma2[k])
# Accept-reject
if (l.metrop + l.hastings > sigma2.thresh[k]) {
previous.theta[[1]][k] <- sigma2.p[k]
}
}
tau = previous.theta[[2]]
# updating tau
tau.p <- exp(log(tau) + tau.jump)
for(k in 1:K)
{
theta = previous.theta
theta[[2]][k] = tau.p[k]
l.metrop <- Log.post.theta(K, od, theta) - Log.post.theta(K, od, previous.theta)
l.metrop <- l.metrop - tau.p[k] * tau.prior.b + (tau.prior.a -1)*log(tau.p[k]) + tau[k] *tau.prior.b - (tau.prior.a -1)*log(tau[k])
l.hastings <- log(tau.p[k]) - log(tau[k])
# Accept-reject
if (l.metrop + l.hastings > tau.thresh[k]) {
previous.theta[[2]][k] <- tau.p[k]
}
}
# updating rho
if(K > 1)
{
rho = previous.theta[[3]]
inv.cdf <- runif(1, min = pnorm(-1, mean = rho, sd = rho.proposal.scale),
max = pnorm(1, mean = rho, sd = rho.proposal.scale))
rho.p <- qnorm(inv.cdf, mean = rho, sd = rho.proposal.scale)
theta = previous.theta
theta[[3]] = rho.p
l.metrop <- Log.post.theta(K, od, theta) - Log.post.theta(K, od, previous.theta)
l.hastings <- -log(pnorm((1 - rho.p) / rho.proposal.scale) - pnorm((-1 - rho.p) / rho.proposal.scale)) +
log(pnorm((1 - rho) / rho.proposal.scale) - pnorm((-1 - rho) / rho.proposal.scale))
# Accept-reject
if (l.metrop + l.hastings > rho.thresh) {
previous.theta[[3]] <- rho.p
}
}
return(previous.theta)
}
#' posterior sampling of beta
postBeta <- function(od, X, K, theta, micro, hyperparam, timesca = 100, check = FALSE) # X: T *K
{
time.d = od$time.comb
filter_ind = od$filter_ind
lc.comb = od$lc + od$c.pred # original observations
if(check) print(lc.comb[c(1,12,23,37,42)])
c.pred = od$c.pred # include mu
rho = theta[[3]]
sigma = sqrt(theta[[1]])
tau = theta[[2]]
leng.time = length(time.d)
for(t in 1:leng.time)
X[t, ceiling(filter_ind[t]/2)] = X[t, ceiling(filter_ind[t]/2)] + c.pred[t]
delta.t <- diff(time.d)
a.i <- exp( -delta.t %*% t(tau)) # (T-1) * K
time.d.t = time.d/timesca
if (micro == 0) {
mat.temp <- matrix(c(rep(1, leng.time)), ncol = 1)
} else if (micro == 1) {
mat.temp <- matrix(c(rep(1, leng.time), time.d.t), ncol = 2)
} else if (micro == 2) {
mat.temp <- matrix(c(rep(1, leng.time), time.d.t, time.d.t^2), ncol = 3)
} else if (micro == 3) {
mat.temp <- matrix(c(rep(1, leng.time), time.d.t, time.d.t^2, time.d.t^3), ncol = 4)
} else if (micro == 4) {
mat.temp <- matrix(c(rep(1, leng.time), time.d.t, time.d.t^2, time.d.t^3, time.d.t^4), ncol = 5)
} else if (micro == 5) {
mat.temp <- matrix(c(rep(1, leng.time), time.d.t, time.d.t^2, time.d.t^3, time.d.t^4, time.d.t^5), ncol = 6)
}
Y = X[-1,] - a.i * X[-leng.time,] # T * K
Y = rbind(X[1,] , Y)
Q = t(rho * sigma) * sigma
tauM = matrix(0, K, K)
for(i in 1:K)
{
for(j in 1:K){
tauM[i, j] = tau[i] + tau[j]
}
}
Q = Q /tauM # Q
micro = micro + 1
Sigma_B = diag(rep(hyperparam[[2]], 2*K)) #original 1/hyperparam[[2]] , take the inverse after
mB = matrix(rep(0, 2* micro * K),ncol=1)
Sigma_B = cbind(Sigma_B, matrix(0, nrow(Sigma_B), K-1))
Sigma_B = rbind(Sigma_B, matrix(0, K-1, ncol(Sigma_B)))
Sigma_B[(2*micro*K + 1) : nrow(Sigma_B), (2*micro*K + 1) : nrow(Sigma_B)] = diag(rep(1, K-1), K-1,K-1)
mB = rbind(mB, matrix(rep(0, K-1), K-1,1))
for(t in 1:leng.time)
{
k = ceiling(filter_ind[t]/2)
f = filter_ind[t]%%2
W = matrix(0, K, 2 * micro * K)
I = matrix(0, K, K - 1)
if(f)
{
W1 = c(mat.temp[t, ] , rep(0,micro))
}else{
W1 = c(rep(0,micro), mat.temp[t, ] ) # 1 * 2m
}
if(t>1)
{
if(delta.t[t-1] ==0) next
k_1 = ceiling(filter_ind[t-1]/2)
f_1 = filter_ind[t-1]%%2
if(f_1)
{
W2 = c(-a.i[t-1, k_1] * mat.temp[t-1, ] , rep(0,micro))
}else{
W2 = c(rep(0,micro), -a.i[t-1, k_1] * mat.temp[t-1, ] ) # 1 * 2m
}
Qt = Q * (1 - exp(-tauM*delta.t[t-1]))
a = a.i[t-1, k_1]
}else{
Qt = Q
W2 = rep(0, 2 * micro)
k_1 = k
f_1 = f
a = 0
}
inv_Qt = chol2inv(chol(Qt))
range = (2 * (k-1) * micro +1) : (2 * k * micro)
range1 = (2 * (k_1-1) * micro +1) : (2 * k_1 * micro)
W[k, range] = W[k, range] + W1
W[k_1, range1] = W[k_1, range1] + W2
if(k > 1) I[k, k-1] = I[k, k-1] + 1
if(k_1 > 1) I[k_1, k_1 -1] = I[k_1, k_1 -1] - a
W = cbind(W, I)
#print(min(eigen(Sigma_B, symmetric = T)$value))
Sigma_B = Sigma_B + t(W)%*%inv_Qt%*%W
mB= mB + t(Y[t, ] %*% inv_Qt %*% W)
}
#print(eigen(Sigma_B, symmetric = T)$value)
Sigma_B = chol2inv(chol(Sigma_B))
mB = Sigma_B %*% mB
B_comb = rmvnorm(1, mB, Sigma_B)
if(any(is.na(B_comb)))
{
B_comb = t(mB)
}
beta = list()
range=0
for(k in 1:K)
{
beta[[k]] = matrix(B_comb[(range + 1) : (range + 2*micro)], ncol=2)
range = range + 2*micro
}
c.pred <- sapply(1:leng.time, function(i) # T * 1 , will change if delta change!!!
{
k = ceiling(filter_ind[i]/2)
f = filter_ind[i] %% 2
if(k>1) {
return(mat.temp[i,] %*% beta[[k]][,2 - f] + B_comb[1, 2*micro*K +k -1]) # T * K = T * m , m * K
}else{
return(mat.temp[i,] %*% beta[[k]][,2 - f])
}
})
lc.comb <- lc.comb - c.pred
mu = c(0, B_comb[1,-1:(-2*micro*K)])
return(list("beta" = beta, "mu" = mu, "lc.comb" = lc.comb, "c.pred" = c.pred))
}
Log.post.delta <- function(delta, data_all, theta, c, c.mu, unif, micro) {
# theta: list, sigma, tau each is a K-vec, rho: K * K matrix
# c: list, beta, each is a m*2 matrix for m-degree poly and 2 curves
K = length(data_all) # number of filters
if (delta < unif[1] | delta > unif[2]) {
-Inf
}
sigma2 <- theta[[1]] # vector
tau <- theta[[2]] # actually is inverse tau!!
rho <- theta[[3]] # K * K matrix: (1, rho, rho,1)
od = orderData(data_all, delta, c, c.mu, micro = micro)
time.comb = od$time.comb
filter_ind = od$filter_ind
var.lc.comb = od$var
lc.comb = od$lc
leng.time = length(time.comb)
# omega.i, i = 1, 2, ..., 2n to be saved
Sigma <- list() # T: K*K, sigma t|t-1
Sigma.star <- list() # T: K*K, sigma t|t
Sigma[[1]] = t(rho * sqrt(sigma2)) * sqrt(sigma2)
temp = matrix(0, K, K)
for(i in 1:K)
{
for(j in 1:K){
temp[i, j] = tau[i] + tau[j]
}
}
Sigma[[1]] = Sigma[[1]] /temp # Q
# x.hat.i, i = 1, 2, ..., 2n to be saved
mu.i <- matrix(0, leng.time, K) # T * K, mu t|t-1
mu.star.i <- matrix(0, leng.time, K) # mu t|t
# a.i, i = 2, ..., 2n
delta.t <- diff(time.comb)
a.i <- exp( -delta.t %*% t(tau)) # (T-1) * K
for (t in 1 : leng.time) {
if(t>1) Sigma[[t]] <- t(a.i[t - 1, ] *Sigma.star[[t-1]]) * a.i[t - 1, ] + Sigma[[1]] * (1 - exp(-temp*delta.t[t-1]))
k = ceiling(filter_ind[t]/2)
Sigma.star[[t]] <- Sigma[[t]] - Sigma[[t]][, k] %*% t(Sigma[[t]][,k]) /(Sigma[[t]][k, k] + var.lc.comb[t])
if(t>1) mu.i[t,] <- a.i[t - 1, ] * mu.star.i[t-1, ]
mu.star.i[t,] <- mu.i[t,] + Sigma[[t]][, k] /(Sigma[[t]][k, k] + var.lc.comb[t]) * (lc.comb[t] - mu.i[t,k])
}
# log-likelihood
mus <- sapply(1:leng.time, function(t) mu.i[t, ceiling(filter_ind[t]/2)])
vars <- sapply(1:leng.time, function(t) Sigma[[t]][ceiling(filter_ind[t]/2), ceiling(filter_ind[t]/2)] + var.lc.comb[t])
if(any(vars<0))
{
print(vars[vars<0])
print(time.comb[which(vars<0)])
}
return(sum(dnorm(lc.comb, mean = mus, sd = sqrt(vars), log = TRUE)))
}
postX.multi <- function(od, K, theta) { # input combined data remove polynomial
# theta: list, sigma2, tau each is a K-vec, rho K * K matrix
time.comb = od$time.comb
filter_ind = od$filter_ind
var.lc.comb = od$var
lc.comb = od$lc
leng.time.comb = length(time.comb)
sigma2 <- theta[[1]] # vector
sigma = sqrt(sigma2)
tau <- theta[[2]] # actually is inverse tau!!
rho <- theta[[3]] # K * K matrix: (1, rho, rho,1)
Sigma <- list() # T: K*K, sigma t|t-1
Sigma.star <- list() # T: K*K, sigma t|t
Sigma[[1]] = t(rho * sigma) * sigma
temp = matrix(0, K, K)
for(i in 1:K)
{
for(j in 1:K){
temp[i, j] = tau[i] + tau[j]
}
}
Sigma[[1]] = Sigma[[1]] /temp # Q
# x.hat.i, i = 1, 2, ..., 2n to be saved
mu.i <- matrix(0, leng.time.comb, K) # T * K, mu t|t-1
mu.star.i <- matrix(0, leng.time.comb, K) # mu t|t
# a.i, i = 2, ..., 2n
delta.t <- diff(time.comb)
a.i <- exp( -delta.t %*% t(tau)) # (T-1) * K
for (t in 1 : leng.time.comb) {
if(t>1) Sigma[[t]] <- t(a.i[t - 1, ] *Sigma.star[[t-1]]) * a.i[t - 1, ] + Sigma[[1]] * (1 - exp(-temp*delta.t[t-1]))
k = ceiling(filter_ind[t]/2)
Sigma.star[[t]] <- Sigma[[t]] - Sigma[[t]][, k] %*% t(Sigma[[t]][,k]) /(Sigma[[t]][k, k] + var.lc.comb[t])
if(t>1) mu.i[t,] <- a.i[t - 1, ] * mu.star.i[t-1, ]
mu.star.i[t,] <- mu.i[t,] + Sigma[[t]][, k] /(Sigma[[t]][k, k] + var.lc.comb[t]) * (lc.comb[t] - mu.i[t,k])
}
# backward sampling
X = matrix(0, leng.time.comb, K)
mut = matrix(0, leng.time.comb, K)
X[leng.time.comb,] <- rmvnorm(1, mu.star.i[leng.time.comb,], Sigma.star[[leng.time.comb]])
mut[leng.time.comb,] <- mu.star.i[leng.time.comb,]
for (t in (leng.time.comb-1) :1) {
mu_t = mu.star.i[t,] + Sigma.star[[t]] %*% ( chol2inv(chol(Sigma[[t+1]])) * a.i[t,])%*% (X[t+1, ] - a.i[t,]* mu.star.i[t,])
sigma_t = Sigma.star[[t]] - Sigma.star[[t]] %*% (chol2inv(chol(Sigma[[t+1]]))* a.i[t,]) %*% (a.i[t,] * Sigma.star[[t]])
mut[t, ] = mu.star.i[t,] + Sigma.star[[t]] %*% ( chol2inv(chol(Sigma[[t+1]])) * a.i[t,])%*% (mut[t+1, ] - a.i[t,]* mu.star.i[t,])
if(delta.t[t] ==0)
{
X[t, ] = X[t+1, ]
}else{
X[t, ]<- rmvnorm(1, mu_t, sigma_t)
}
}
#X + c.pred
return(list('X' = X, 'mu' = mut))
}
#' @param c: microlensing polynomial coefficients
#' @param c.mu: mean diff of different bands
#' @return c.pred: polynomial part of the observations(including the means)
#' @return lc.comb: residual part
orderData <- function(data_all, delta, c, c.mu, timesc = 100, remove.poly = TRUE, micro = NULL)
{
# sorting time given delta
time.temp <- c()
filter_ind <- c()
var.lc.comb <- c()
lc.comb <- c() # observations
k = 1
for(data in data_all)
{
time <- data[, 1]
time.d <- time - delta
time.temp <- c(time.temp, time, time.d)
filter_ind <-c(filter_ind, rep(k, length(time)), rep(k+1, length(time)))
var.lc.comb <- c(var.lc.comb , data[,3]^2, data[,5]^2)
lc.comb <- c(lc.comb, data[,2], data[,4])
k = k + 2
}
ord <- order(time.temp)
time.comb <- time.temp[ord]
filter_ind <- filter_ind[ord]
var.lc.comb <- var.lc.comb[ord]
lc.comb <- lc.comb[ord]
if(!remove.poly) return(list("time.comb" = time.comb, "filter_ind" = filter_ind, "var" = var.lc.comb, "lc" = lc.comb))
leng.time <- length(time.comb)
time.comb.t = time.comb/timesc
if (micro == 0) {
mat.temp <- matrix(c(rep(1, leng.time)), ncol = 1)
} else if (micro == 1) {
mat.temp <- matrix(c(rep(1, leng.time), time.comb.t), ncol = 2)
} else if (micro == 2) {
mat.temp <- matrix(c(rep(1, leng.time), time.comb.t, time.comb.t^2), ncol = 3)
} else if (micro == 3) {
mat.temp <- matrix(c(rep(1, leng.time), time.comb.t, time.comb.t^2, time.comb.t^3), ncol = 4)
}
c.pred <- sapply(1:leng.time, function(i) # T * 1 , will change if delta change!!!
{
k = ceiling(filter_ind[i]/2)
f = filter_ind[i] %% 2
return(mat.temp[i,] %*% c[[k]][,2 - f] + c.mu[k]) # T * K = T * m , m * K
})
lc.comb <- lc.comb - c.pred
return(list("time.comb" = time.comb, "filter_ind" = filter_ind, "var" = var.lc.comb, "lc" = lc.comb, "c.pred" = c.pred))
}
#' mcmc sampling
#' @param data_all: a list, each element is a data frame for one band, columns are time, value and std for different light curves (delayed).
#' @param theta.ini: initial values of the OU process parameters of multiple bands
#' @param delta.ini: initial value for time delay
###' @param kappa.ini: initial hyperparameter for the prior of microlensing polynomial coefficients
#' @param c.ini: initial coefficients of polynomials. A list, each element is a data frame for one band. Rows are coefficients of different degrees, columns are different light curves. Default: NULL
#' @param delta.uniform.range: prior of possible range of time delay
#' @param tau.*: paramters for prior, hyperprior and mcmc proposal distribution of 1/tau in OU process
#' @param sigma.*: paramters for prior, hyperprior and mcmc proposal distribution of sigma in OU process
#' @param rho.*: paramters for prior and mcmc proposal distribution of rho in OU process
###' @param kappa.*: paramters for prior and mcmc proposal distribution of kappa in the prior of microlensing polynomial coefficients
###' @param hyper.sigma2.*: hyperparameters for the microlensing polynomial coefficients
#' @param beta.prior.diag: prior for microlensing polynomial coefficients
#' @param micro: the degree of polynomial for microlensing
#' @param g: scale the time for fitting polynomial of microlensing
#' @param adaptive.*: parameters for adaptive MCMC
#' @param sample.size, warmingup,size: mcmc sample and burnin size
bayesian.multiband <- function(data.band1, data.band2, n.bands = 2,
theta.ini = c(0.01, 0.01, 100, 100, 0.5),
delta.ini, delta.uniform.range = c(-500, 500), delta.proposal.scale = 1,
tau.proposal.scale = 1, tau.prior.shape = 1, tau.prior.scale = 1,
sigma2.proposal.scale = 0.5, sigma2.prior.shape = 1, sigma2.prior.scale = 1e-7,
rho.proposal.scale = 0.1,
beta.prior.diag = 10 * c(0.1, 0.01, 1e-3, 1e-5)^2,
micro = 3, timesc = 100, adaptive.frequency = 100,
adaptive.delta.factor = 0.1, adaptive.tau.factor = 0.1,
adaptive.sigma2.factor = 0.1, adaptive.rho.factor = 0.1,
sample.size, warmingup.size) {
min.time <- min(data.band1[1, 1], data.band2[1, 1])
data.band1[, 1] <- data.band1[, 1] - 3500
data.band2[, 1] <- data.band2[, 1] - 3500
K <- n.bands
data_all <- vector(mode = "list", length = K)
for (i in 1 : K) {
eval(parse(text = paste0("data.temp <- data.band", i)))
colnames(data.temp)[1] <- "V1"
data_all[[i]] <- data.temp
}
beta.prior.diag <- beta.prior.diag[1 : (micro + 1)]
total.sample.size <- sample.size + warmingup.size
leng.time.comb = sum(sapply(data_all, nrow))
sigma2.out <- matrix(NA, total.sample.size, K)
tau.out <- matrix(NA, total.sample.size, K)
hyper.prior.out <- matrix(NA, total.sample.size, 4) # tau.shape, tau.scale, sigma2.shape, sigma2.scale
delta.out <- rep(NA, total.sample.size)
rho.out <- rep(NA, total.sample.size)
c.out <- matrix(NA, nrow = total.sample.size, ncol = (micro + 1) * 2 * K)
c.mu.out <- matrix(NA, nrow = total.sample.size, ncol = K)
X.out <- matrix(NA, nrow = total.sample.size, ncol = leng.time.comb * K*2)
loglik.out <- rep(NA, total.sample.size)
tau.accept <- matrix(0, total.sample.size, K)
sigma2.accept <- matrix(0, total.sample.size, K)
delta.accept <- rep(0, total.sample.size)
rho.accept <- rep(0, total.sample.size)
tau.jumps <- tau.proposal.scale * matrix(rnorm(total.sample.size * K), nrow = K, ncol =total.sample.size)
tau.thresh <- matrix(-rexp(total.sample.size * K), ncol = K, nrow =total.sample.size)
tau.proposal.scale.adapt <- rep(1,K)
sigma2.jumps <- sigma2.proposal.scale * matrix(rnorm(total.sample.size * K), nrow = K, ncol =total.sample.size)
sigma2.thresh <- matrix(-rexp(total.sample.size * K), ncol = K, nrow =total.sample.size)
sigma2.proposal.scale.adapt <- rep(1,K)
delta.jumps <- delta.proposal.scale * rnorm(total.sample.size)
delta.thresh <- -rexp(total.sample.size)
delta.proposal.scale.adapt <- 1
rho.thresh <- -rexp(total.sample.size)
rho.proposal.scale.adapt <- rho.proposal.scale
delta.t <- delta.ini # delta ini
sigma2.t <- theta.ini[1 : 2]^2 #sigma ini
tau.t <- 1 / theta.ini[3 : 4] #tau ini
rho <- theta.ini[5]
if(K > 1){
rho.t = matrix(rho, nrow = K, ncol =K)
diag(rho.t) = 1
}else rho.t = matrix(rho)
theta.prior.t = list(tau.prior.shape, tau.prior.scale, sigma2.prior.shape, sigma2.prior.scale)
hyperparam.t = list(0, beta.prior.diag, rep(0, micro+1))
c.t = list()
c.mu.t = rep(0, K)
i =1
mu0 = mean(c(data_all[[1]][,2], data_all[[1]][,4])) # mean of different light curves for the first band
for(dat in data_all)
{
c.mu.t[i] = mean(c(dat[,2], dat[,4])) - mu0 # different in mean relatively to the first band
ti1 <- dat[, 1]/timesc
ti2 <- ti1^2
ti3 <- ti1^3
if (micro == 0) {
ss1 <- lm(dat[, 2] - c.mu.t[i] ~ 1)
} else if (micro == 1) {
ss1 <- lm(dat[, 2] - c.mu.t[i] ~ ti1)
} else if (micro == 2) {
ss1 <- lm(dat[, 2] - c.mu.t[i] ~ ti1 + ti2)
} else if (micro == 3) {
ss1 <- lm(dat[, 2] - c.mu.t[i] ~ ti1 + ti2 + ti3)
}
ti1 <- (dat[, 1] - delta.ini)/timesc
ti2 <- ti1^2
ti3 <- ti1^3
if (micro == 0) {
ss2 <- lm(dat[, 4] - c.mu.t[i] ~ 1)
} else if (micro == 1) {
ss2 <- lm(dat[, 4] - c.mu.t[i] ~ ti1)
} else if (micro == 2) {
ss2 <- lm(dat[, 4] - c.mu.t[i] ~ ti1 + ti2)
} else if (micro == 3) {
ss2 <- lm(dat[, 4] - c.mu.t[i] ~ ti1 + ti2 + ti3)
}
if(micro > 3){
c.t[[i]] <- cbind(c(ss1$coefficients, rep(0,micro - 2)), c(ss2$coefficients, rep(0,micro - 2)))
}else{
c.t[[i]] <- cbind(ss1$coefficients, ss2$coefficients)
}
i = i +1
}
print(paste("Starting time:", Sys.time()))
od <- orderData(data_all, delta.t, c.t, c.mu.t, remove.poly = T, micro = micro) # remove the means of different bands/curves
for (i in 1 : total.sample.size) {
reslist <- postX.multi(od, K, theta = list(sigma2.t, tau.t, rho.t))
X.t = reslist$X
mu.t = reslist$mu
tmp_list = postBeta(od, X.t, K, theta = list(sigma2.t, tau.t, rho.t), micro, hyperparam.t, timesca =timesc, check = FALSE) # X: T *K
c.t = tmp_list$beta
od$lc.comb = tmp_list$lc.comb
c.mu.t = tmp_list$mu
# theta update
tau.jump.adapt <- tau.proposal.scale.adapt * tau.jumps[,i]
sigma2.jump.adapt <- sigma2.proposal.scale.adapt * sigma2.jumps[,i]
if(K > 1) rho = rho.t[1,2]
theta.update <- postTheta.multi(K, od, previous.theta = list(sigma2.t, tau.t, rho),
tau.jump = tau.jump.adapt, tau.thresh = tau.thresh[i, ],
tau.prior.a = tau.prior.shape, tau.prior.b = tau.prior.scale,
sigma2.jump = sigma2.jump.adapt, sigma2.thresh = sigma2.thresh[i, ],
sigma.prior.a = sigma2.prior.shape, sigma.prior.b = sigma2.prior.scale,
rho.proposal.scale = rho.proposal.scale.adapt, rho.thresh = rho.thresh[i])
if (K > 1){
if( theta.update[[3]] != rho.t[1,2]) {
rho.accept[i] <- 1
rho.t = matrix(theta.update[[3]], nrow = K, ncol =K)
diag(rho.t) = 1
}
}
tau.accept[i, theta.update[[2]] != tau.t] <- 1
sigma2.accept[i, theta.update[[1]] != sigma2.t] <- 1
# delta update
delta.p <- delta.t + delta.proposal.scale.adapt * delta.jumps[i]
prev_lik <- Log.post.delta(delta.t, data_all, list(sigma2.t, tau.t, rho.t), c.t, c.mu.t,
unif = delta.uniform.range, micro)
cur_lik <- Log.post.delta(delta.p, data_all, list(sigma2.t, tau.t, rho.t), c.t, c.mu.t,
unif = delta.uniform.range, micro)
l.metrop <- cur_lik - prev_lik
if (l.metrop > delta.thresh[i]) {
delta.t <- delta.p
delta.accept[i] <- 1
od <- orderData(data_all, delta.t, c.t, c.mu.t, remove.poly = TRUE, micro = micro)
loglik.out[i] <- cur_lik
}else loglik.out[i] <- prev_lik
sigma2.t <- sigma2.out[i,] <- theta.update[[1]]
tau.t <- tau.out[i,] <- theta.update[[2]]
rho.out[i] <- theta.update[[3]]
X.out[i,] <- mu.t
c.out[i, ]<- unlist(c.t)
c.mu.out[i, ] <- c.mu.t
delta.out[i] <- delta.t
if (i %% adaptive.frequency == 0) {
if(mean(delta.accept[i - (adaptive.frequency - 1) : i]) > 0.35) {
scale.adj <- exp(min(adaptive.delta.factor, 1 / sqrt(i / adaptive.frequency)))
} else if (mean(delta.accept[i - (adaptive.frequency - 1) : i]) < 0.35) {
scale.adj <- exp(-min(adaptive.delta.factor, 1 / sqrt(i / adaptive.frequency)))
} else {
scale.adj <- 1
}
delta.proposal.scale.adapt <- delta.proposal.scale.adapt * scale.adj
}
if (i %% adaptive.frequency == 0) {
if(mean(rho.accept[i - (adaptive.frequency - 1) : i]) > 0.35) {
scale.adj <- exp(min(adaptive.delta.factor, 1 / sqrt(i / adaptive.frequency)))
} else if (mean(rho.accept[i - (adaptive.frequency - 1) : i]) < 0.35) {
scale.adj <- exp(-min(adaptive.delta.factor, 1 / sqrt(i / adaptive.frequency)))
} else {
scale.adj <- 1
}
rho.proposal.scale.adapt <- rho.proposal.scale.adapt * scale.adj
}
if (i %% adaptive.frequency == 0) {
scale.adj <- exp(min(adaptive.tau.factor, 1 / sqrt(i / adaptive.frequency)))
ind = which(colMeans(tau.accept[(i - (adaptive.frequency - 1)) : i, , drop=FALSE]) > 0.35)
tau.proposal.scale.adapt[ind] <- tau.proposal.scale.adapt[ind] * scale.adj
ind = which(colMeans(tau.accept[(i - (adaptive.frequency - 1)) : i, , drop=FALSE]) < 0.35)
tau.proposal.scale.adapt[ind] <- tau.proposal.scale.adapt[ind] / scale.adj
}
if (i %% adaptive.frequency == 0) {
scale.adj <- exp(min(adaptive.sigma2.factor, 1 / sqrt(i / adaptive.frequency)))
ind = which(colMeans(sigma2.accept[(i - (adaptive.frequency - 1)) : i, , drop=FALSE]) > 0.35)
sigma2.proposal.scale.adapt[ind] <- sigma2.proposal.scale.adapt[ind] * scale.adj
ind = which(colMeans(sigma2.accept[(i - (adaptive.frequency - 1)) : i, , drop=FALSE]) < 0.35)
sigma2.proposal.scale.adapt[ind] <- sigma2.proposal.scale.adapt[ind] / scale.adj
}
}
print(paste("Ending time:", Sys.time()))
out <- list(delta = delta.out[-c(1 : warmingup.size)],
beta = c.out[-c(1 : warmingup.size), ],
rho = rho.out[-c(1 : warmingup.size)],
sigma = sqrt(sigma2.out[-c(1 : warmingup.size),]),
tau = 1 / tau.out[-c(1 : warmingup.size),],
tau.accept.rate = colMeans(tau.accept),
sigma.accept.rate = colMeans(sigma2.accept),
delta.accept.rate = mean(delta.accept),
rho.accept.rate = mean(rho.accept))
out
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.