R/rmixmorm.R

In thiyangt/fformpp: Feature-based FORecast Model Performance Prediction

##' Generate random variables from univariate and multivariate mixture normal distribution
##'
##' @export
rmixnorm <- function(n, means, sigmas, weights)
{

    if(!is.matrix(means) | length(dim(sigmas)) != 3)
    {
        stop("means must be a q-by-k matrix and sigmas must be a q-by-q-by-k array.")
    }

    k <- length(weights) # K-components
    q <- dim(means)[1] # q-dimensional

    idx <- rmultinom(n = n, 1, prob = weights) # k-by-n matrix

    out <- apply(idx, 2, function(x, means, sigmas, q)
    {
        which.comp <- which(x == 1)
        mvtnorm::rmvnorm(n = 1, mean = means[, which.comp],
                sigma = matrix(sigmas[, , which.comp], q, q))
    }, means = means, sigmas = sigmas, q = q)

    return(out)
}



##' Density for univariate and multivariate mixture normal distribution
##'
##' Density from mixture of normals.
##' @title Density for normal mixtures
##' @param n An positive integer. Numbers of samples to be generated.
##' @param means A `q-by-k` matrix. Mean value vector within each component in total k
##'     components.
##' @param sigmas A `q-by-q-by-k` array.  Variance covariance matrix with in each
##'     component.
##' @param weights A `k`-length vector.  The weight in each component.
##' @return An `n-by-q` matrix
##' @references Villani et al 2009
##' @author Feng Li, Central University of Finance and Economics.
##' @export
dmixnorm <- function(x, means, sigmas, weights, log = FALSE)
{

    if(!is.matrix(means) | length(dim(sigmas)) != 3)
    {
        stop("means must be a q-by-k matrix and sigmas must be a q-by-q-by-k array.")
    }

    k <- length(weights) # K-components
    q <- dim(means)[1] # q-dimensional

    out.comp.log <- apply(matrix(1:k), 1, function(comp.i, x, means, sigmas, q)
    {
        mvtnorm::dmvnorm(x = matrix(x, 1, q), mean = matrix(means[, comp.i], 1, q),
                sigma = matrix(sigmas[, ,comp.i], q, q),
                log = TRUE)
    }, x = x,  means = means, sigmas = sigmas, q = q)


    out.sum <- sum(exp(out.comp.log)*weights)

    if(log  == TRUE)
    {
        out <- log(out.sum)
    }
    else
    {
        out <- out.sum
    }

    return(out)
}



##' Simulating AR type random variables from mixture of normals
##'
##' This function simulates random samples from a finite mixture of Gaussian distribution
##'     where the mean from each components are AR(p) process.
##' @title AR-type random variables from mixture of normals
##' @param n a positive integer.  Number of samples.
##' @param means.ar.par.list parameters in AR(p) within each mixing component
##' @param sigmas.list variance list
##' @param weights weight in each list
##' @param yinit initial values
##' @return vector of n follows a mixture distribution
##' @references Li 2010 JSPI
##' @author Feng Li, Central University of Finance and Economics.
##' @export
rmixnorm.ts <- function(n, means.ar.par.list, sigmas.list, weights, yinit = 0)
{
    y <- rep(NA, n)
    nComp <- length(means.ar.par.list)
    nLags <- lapply(means.ar.par.list, function(x) length(x)-1)
    maxLag <- max(unlist(nLags))

    if(any(unlist(nLags)<1))
    {
        stop("Drift is always included. Set the first elements in ar.par.list be zero to remove the drift effect.")
    }

    y[1:maxLag] <- yinit

    sigmas.ary <- matrix(do.call(cbind, sigmas.list), n, nComp, byrow = TRUE)

    for(i in (maxLag+1):n)
    {
        meansComp <- lapply(means.ar.par.list, function(par, y, i){
            nLag <- length(par)-1
            yPre <- c(1, y[(i-1):(i-nLag)])
            yCurr <- sum(par*yPre)
            ## print(cbind(yPre, par))
            return(yCurr)
        }, y = y, i = i)
        sigmas.i <- array(sigmas.ary[i, ], dim = c(1, 1, nComp))
        y[i] <- mvtnorm::rmixnorm(n = 1, means = matrix(unlist(meansComp), nrow = 1),
                         sigmas = sigmas.i, weights = weights)
    }
    return(as.ts(y))
}


##' @title The conditional density for AR-type mixture of normals
##' @param y An `n`-length vector. The time series.
##' @param means.ar.par.list An `k`-length list. Each element in the list consists the
##'     coefficients in the AR process.
##' @param sigmas.list
##' @param weights A `k`-length vector. The weight in each component.
##' @param log Logical; If TRUE, the log density is returned.
##' @export
dmixnorm.ts <- function(y, means.ar.par.list, sigmas.list, weights, log = FALSE)
{
    nComp <- length(means.ar.par.list)
    nLags <- lapply(means.ar.par.list, function(x) length(x)-1)
    maxLag <- max(unlist(nLags))
    n <- length(y)

    if(any(unlist(nLags)<1))
    {
        stop("Drift is always included. Set the first elements in ar.par.list be zero to remove the drift effect.")
    }

    out.log <- y
    out.log[] <- 0 # Set the density be zero for the first maxlags.

    sigmas.ary <- matrix(do.call(cbind, sigmas.list), n, nComp, byrow = TRUE)

    for(i in (maxLag+1):n)
    {
        meansComp <- lapply(means.ar.par.list, function(par, y, i){
            nLag <- length(par)-1
            yPre <- c(1, y[(i-1):(i-nLag)])
            yCurr <- sum(par*yPre)
            ## print(cbind(yPre, par))
            return(yCurr)
        }, y = y, i = i)
        sigmas.i <- array(sigmas.ary[i, ], dim = c(1, 1, nComp))
        out.log[i] <- rmvtnorm::dmixnorm(x = y[i], means = matrix(unlist(meansComp), nrow = 1),
                               sigmas = sigmas.i, weights = weights, log = TRUE)
    }

    if(log == TRUE)
    {
        out <- out.log
    }
    else
    {
        out <- exp(log)
    }
    return(out)
}