R/marg_loglik.R

In mixAR: Mixture Autoregressive Models

## 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)
}