## Do not edit this file manually.
## It has been automatically generated from mixAR.org.
marg_loglik <- function(y, model, tau, nsim, prob_mod){
## Note: "prob_mod" is the posterior probability of choosing the model in the
## RJMCMC procedure ("choosepk")-- i.e. the largest output of the function
## should approximately give p(p*|y)
if(prob_mod > 1 | prob_mod < 0) stop("prob_mod should be between 0 and 1")
n <- length(y)
pk <- model@order
p <- max(pk)
g <- length(pk)
AR <- model@arcoef@a
bk <- sapply(AR, function(x) 1 - sum(x))
mu <- model@shift;
mu[which(model@shift != 0)] <- model@shift[which(model@shift != 0)] / bk[which(model@shift != 0)]
rat <- matrix(nrow=nsim, ncol=g)
if(length(tau) == 1) tau <- rep(tau, g)
if(length(tau) != 1 & length(tau) != g){
stop("tau should have length equal to 1 or g")
}
lnames <- NULL
ard <- list()
phi_hd <- list()
for(k in 1:g){
ard[[k]] <- matrix(0,nrow=nsim,ncol=pk[k])
phi_hd[[k]] <- numeric(pk[k])
lnames <- c(lnames, paste("Component ", k))
}
names(ard) <- lnames
## Hyperparameters
d <- rep(1,g)
a <- .2
c <- 2
Ry <- max(y) - min(y)
zeta <- min(y) + Ry / 2
ka <- 1 / Ry
b <- 100 * a / (c * Ry^2)
## Begin simulation to estimate phi_hd
## The remaining parameters here are free to vary
count <- rep(0,g)
for(j in 1:nsim){
if(j == 1){
xx <- sampZpi(y, pk, model@prob, mu, AR, model@scale, 1, d)
xx2 <- sampMuShift(y, pk, 1/model@scale^2, xx$nk, model@shift,
xx$latentZ, AR, 1)
xx3 <- sampSigmaTau(y, pk, 1/model@scale^2, xx$nk, AR, xx2$mu,
xx$latentZ, a, c, 1)
}else{
xx <- sampZpi(y, pk, xx$mix_weights, xx2$mu, AR, xx3$scale, 1, d)
xx2 <- sampMuShift(y, pk, xx3$precision, xx$nk, xx2$shift,
xx$latentZ, AR, 1)
xx3 <- sampSigmaTau(y, pk, xx3$precision, xx$nk, AR, xx2$mu,
xx$latentZ, a, c, 1)
}
ratio <- rep(0, g)
phi_k <- list()
for(k in 1:g){
phi_k[[k]] <- rnorm(pk[k], AR[[k]], tau[k])
ratio[k] <- min(1,exp(- xx3$precision[k]/2 *
sum(err_k(phi_k[[k]], xx2$mu[k], y, xx$latentZ[ ,k], p, pk[k]) ^ 2 -
err_k(AR[[k]], xx2$mu[k], y, xx$latentZ[ ,k], p, pk[k]) ^ 2)))
}
cc <- rep(0,g)
u <- runif(g)
rat[j, ] <- ratio
for(k in 1:g){
if(u[k] > ratio[k]){
phi_k[[k]] <- AR[[k]]
}
}
A <- new("MixARGaussian", prob = as.numeric(xx$mix_weights),
scale = as.numeric(xx3$scale), arcoef = phi_k,
shift = as.numeric(xx2$shift))
if(isStable(A)){
for(k in 1:g){
ard[[k]][j, ] <- phi_k[[k]]
}
AR <- phi_k
} else{
for(k in 1:g){
ard[[k]][j, ] <- AR[[k]]
rat[j, k] <- 0
}
}
}
## Set phi1 to its maximum density
for(l in 1:pk[1]){
phi_hd[[1]][l] <- density(ard[[1]][,l])$x[which.max(density(ard[[1]][,l])$y)]
}
## Calculate numerator for estimated density
rr <- 0
for(j in 1:nsim){
ARR <- list(ard[[1]][j,])
rr[j] <- min(1,exp(- xx3$precision[1]/2 *
sum(err_k(phi_hd[[1]], xx2$mu[1], y, xx$latentZ[ ,1], p, pk[1]) ^ 2 -
err_k(ARR[[1]], xx2$mu[1], y, xx$latentZ[ ,1], p, pk[1]) ^ 2)))
}
num <- 0
for(j in 1:nsim){
num[j] <- rr[j] * prod(dnorm(phi_hd[[1]], ard[[1]][j, ], tau[1]))
}
num <- mean(num)
dens_phi <- rep(1, g)
AR[[1]] <- phi_hd[[1]]
## Repeat procedure for k=2,...,g
## Need different routine for k=g, hence the "if"
for(kk in 2:g){
if(kk < g){
ard1 <- list()
for(l in 1:(g - kk + 1)){
ard1[[l]] <- matrix(nrow = nsim, ncol = pk[kk - 1 + l])
}
ratio1 <- numeric(g - kk + 1)
for(j in 1:nsim){
xx <- sampZpi(y, pk, xx$mix_weights, xx2$mu, AR, xx3$scale, 1, d)
xx2 <- sampMuShift(y, pk, xx3$precision, xx$nk, xx2$shift,
xx$latentZ, AR, 1)
xx3 <- sampSigmaTau(y, pk, xx3$precision, xx$nk, AR, xx2$mu,
xx$latentZ, a, c, 1)
phi_tild <- rnorm(pk[kk - 1], phi_hd[[kk - 1]], tau[kk - 1])
## "rr" is denominator for density of component kk-1
rr[j] <- min(1,exp(- xx3$precision[kk-1]/2 *
sum(err_k(phi_tild,
xx2$mu[kk-1], y, xx$latentZ[ ,kk-1], p, pk[kk-1]) ^ 2 -
err_k(phi_hd[[kk-1]],
xx2$mu[kk-1], y, xx$latentZ[ ,kk-1], p, pk[kk-1]) ^ 2)))
phi_k <- phi_hd
for(k in kk:g){
phi_k[[k]] <- rnorm(pk[k], AR[[k]], tau[k])
ratio1[k - kk + 1] <-
min(1, exp(- xx3$precision[k]/2 *
sum(err_k(phi_k[[k]],
xx2$mu[k], y, xx$latentZ[ ,k], p, pk[k])^2 -
err_k(AR[[k]],
xx2$mu[k], y, xx$latentZ[ ,k], p, pk[k])^2)))
}
u <- runif(g - kk + 1)
for(k in 1:(g - kk + 1)){
if(u[k] > ratio1[k]) phi_k[[k + 1]] <- AR[[k + 1]]
}
A <- new("MixARGaussian", prob = as.numeric(xx$mix_weights),
scale = as.numeric(xx3$scale), arcoef = phi_k,
shift = as.numeric(xx2$shift))
if(isStable(A)){
for(l in 1:(g - kk + 1)){
ard1[[l]][j, ] <- phi_k[[l - 1 + kk]] ##l+1
}
AR <- phi_k
}else{
for(l in 1:(g - kk + 1)){
ard1[[l]][j, ] <- AR[[l - 1 + kk]] ##l+1
}
}
}
dens_phi[[kk - 1]] <- num / mean(rr) ## estimate density of the previous component
rr <- 0
num <- 0
## calculate numerator of density for the current component
for(l in 1:pk[kk]){ ##create phi_hd for current component
phi_hd[[kk]][l] <-
density(ard1[[1]][,l])$x[which.max(density(ard1[[1]][,l])$y)]
}
for(j in 1:nsim){
ARR <- list(ard1[[1]][j,])
rr[j] <- min(1,exp(- xx3$precision[kk]/2 *
sum(err_k(phi_hd[[kk]], xx2$mu[kk], y,
xx$latentZ[ ,kk], p, pk[kk]) ^ 2 -
err_k(ARR[[1]], xx2$mu[kk], y, xx$latentZ[ ,kk],
p, pk[kk]) ^ 2)))
num[j] <- rr[j] * prod(dnorm(phi_hd[[kk]], ard1[[1]][j, ], tau[kk]))
}
num <- mean(num)
AR[[kk]] <- phi_hd[[kk]]
}
if(kk == g){
ard1 <- matrix(nrow = (nsim),ncol = pk[g])
rr <- numeric(nsim)
rat1 <- numeric(nsim)
for(j in 1:nsim){
xx <- sampZpi(y, pk, xx$mix_weights, xx2$mu, AR, xx3$scale, 1, d)
xx2 <- sampMuShift(y, pk, xx3$precision, xx$nk, xx2$shift,
xx$latentZ, AR, 1)
xx3 <- sampSigmaTau(y, pk, xx3$precision, xx$nk, AR, xx2$mu,
xx$latentZ, a, c, 1)
phi_tild <- rnorm(pk[kk-1], phi_hd[[kk-1]], tau[kk-1])
rr[j] <- min(1,exp(- xx3$precision[kk-1]/2 *
sum(err_k(phi_tild, xx2$mu[kk-1], y,
xx$latentZ[ ,kk-1], p, pk[kk-1]) ^ 2 -
err_k(phi_hd[[kk-1]], xx2$mu[kk-1], y,
xx$latentZ[ ,kk-1], p, pk[kk-1])^2)))
phi_k <- phi_hd
phi_k[[g]] <- rnorm(pk[g], AR[[g]], tau[g])
ratio1 <- min(1,exp(- xx3$precision[g]/2 *
sum(err_k(phi_k[[g]],
xx2$mu[g], y, xx$latentZ[ ,g], p, pk[g]) ^ 2 -
err_k(AR[[g]],
xx2$mu[g], y, xx$latentZ[ ,g], p, pk[g]) ^ 2)))
u <- runif(1)
rat1[j] <- ratio1
for(k in 1:g){
if(u > ratio1){
phi_k[[g]] <- AR[[g]]
}
A <- new("MixARGaussian", prob = as.numeric(xx$mix_weights),
scale = as.numeric(xx3$scale), arcoef = phi_k,
shift = as.numeric(xx2$shift))
if(isStable(A)){
ard1[j, ] <- phi_k[[g]]
AR <- phi_k
} else {
ard1[j, ] <- AR[[g]]
rat1[j] <- 0
}
}
}
dens_phi[kk - 1] <- num / mean(rr)
phi_hd[[g]] <- numeric(pk[g])
for(l in 1:pk[g]){
phi_hd[[g]][l] <- density(ard1[,l])$x[which.max(density(ard1[,l])$y)]
dens_phi[g] <- dens_phi[g] * max(density(ard1[,l])$y)
}
}
}
## #######################################
## Calculate mu_hd
## phi_hd fixed, remaining parameters free to vary
mud <- matrix(nrow = nsim, ncol = g)
bk <- 0
for(k in 1:g){
bk[k] <- 1- sum(phi_hd[[k]])
}
mu_hd <- numeric(g)
prec <- nk <- dens_mu_m <- matrix(nrow = nsim, ncol = g)
latZ <- list()
for(j in 1:nsim){
xx2 <- sampMuShift(y, pk, xx3$precision, xx$nk, xx2$shift, xx$latentZ,
phi_hd, 1)
mud[j, ] <- xx2$mu
xx <- sampZpi(y, pk, xx$mix_weights, xx2$mu, phi_hd, xx3$scale, 1, d)
nk[j, ] <- xx$nk
latZ[[j]] <- xx$latentZ
xx3 <- sampSigmaTau(y, pk, xx3$precision, xx$nk, phi_hd, xx2$mu,
xx$latentZ, a, c, 1)
prec[j, ] <- xx3$precision
}
for(k in 1:g){
mu_hd[k] <- density(mud[,k])$x[which.max(density(mud[,k])$y)]
for(j in 1:nsim){
dens_mu_m[j, k] <- dnorm(mu_hd[k],(prec[j, k] * nk[j, k] *
mean(err(phi_hd, mud[j, ], y,
latZ[[j]],pk)[, k]) *
bk[k] + ka * zeta)/
(prec[j, k] * nk[j, k] * bk[k]^2 + ka),
1 / sqrt((prec[j, k] * nk[j, k] * bk[k]^2 + ka)))
}
}
dens_mu <- mean(apply(dens_mu_m, 1 ,prod))
bk <- sapply(phi_hd, function(x) 1 - sum(x))
## Calculate shift_hd using formula
shift_hd <- mu_hd * (bk)
## ############################
## Calculation of sigma_hd/prec_hd
## phi_hd & mu_hd fixed, pi_hd variable
precd <- matrix(ncol = g,nrow = nsim)
lambda <- numeric(nsim)
prec_hd <- numeric(g)
dens_prec_m <- matrix(nrow = nsim, ncol = g)
nk <- matrix(nrow = nsim, ncol = g)
for(j in 1:nsim){
xx <- sampZpi(y, pk, xx$mix_weights, mu_hd, phi_hd, xx3$scale, 1, d)
nk[j, ] <- xx$nk
xx3 <- sampSigmaTau(y, pk, xx3$precision, xx$nk, phi_hd, mu_hd,
xx$latentZ, a, c, 1)
precd[j, ] <- xx3$precision
lambda[j] <- xx3$lambda
}
for(k in 1:g){
prec_hd[k] <- density(precd[,k])$x[which.max(density(precd[,k])$y)]
for(j in 1:nsim){
dens_prec_m[j,k] <- dgamma(prec_hd[k],c + nk[j,k]/2, lambda[j] + 0.5 *
sum(err(phi_hd, mu_hd, y, latZ[[j]],
pk)[ ,k]^2))
}
}
## Calculate scale_hd using formula
scale_hd <- 1/ sqrt(prec_hd)
dens_prec <- mean(apply(dens_prec_m,1,prod))
## ###########################
## Calculate pi_hd - all remaining parameters fixed
prob_hd <- numeric(g)
nk <- prob <- matrix(nrow = nsim, ncol = g)
for(j in 1:nsim){
xx <- sampZpi(y, pk, xx$mix_weights, mu_hd, phi_hd, scale_hd, 1, d)
nk[j, ] <- xx$nk
prob[j, ] <- round(xx$mix_weights,3)
}
for(k in 1:(g-1)){
prob_hd[k] <- density(prob[,k])$x[which.max(density(prob[,k])$y)]
}
prob_hd[g] <- 1-sum(prob_hd[1:(g-1)])
dens_prob <- numeric(nsim)
for(j in 1:nsim){
dens_prob[j] <- ddirichlet(prob_hd, d + nk[j, ])
}
## #####
## Calculate marginal likelihood for HD parameters
## For comparison we use loglik..easier to see differences
lik <- matrix(nrow = n-p,ncol = g)
for(t in 1:(n-p)){
for(k in 1:g){
lik[t,k] <- prob_hd[k] * dnorm(y[t+p], mu_hd[k] +
(y[(t + p - 1) : (t+p-pk[k])]-mu_hd[k])%*%
phi_hd[[k]], scale_hd[k])
}
}
lik <- sum(log(rowSums(lik)))
lambda_hd <- density(lambda)$x[which.max(density(lambda)$y)]
priors <- prod(dgamma(prec_hd, c, lambda_hd)) * dgamma(lambda_hd, a, b) *
prod(dnorm(mu_hd,zeta, 1/sqrt(ka)))
for(k in 1:g){
priors <- priors * prod(dnorm(phi_hd[[k]]))
}
posteriors <- prod(dens_phi) * prod(dens_mu) * prod(dens_prec) * mean(dens_prob)
margll <- lik + log(priors) - log(posteriors) + prob_mod
list(marg_loglik = margll, phi_hd = phi_hd, prec_hd = prec_hd,
mu_hd = mu_hd, weig_hd = prob_hd)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.