Computing Mixing Coefficients (Weights) of Mixed Copula Model (Gumbel and Unstructured Student-t) using EM Algorithm in R

Question

Hello Stack Overflow community,

I am currently working on a project involving the estimation of a mixed copula model composed of a Gumbel copula and an unstructured Student-t copula. I would like to know how to compute the mixing coefficients (weights) for these copulas using the Expectation-Maximization (EM) algorithm in R.

More specifically, How do I implement the EM algorithm in R to estimate the mixing coefficients (weights) for Copula Models? I saw that mixtools package should provide some solutions but haven't found the right one,

I have already searched for relevant resources and documentation, but I couldn't find a clear example or tutorial specifically addressing this mixed copula model with the EM algorithm in R. Any guidance, code snippets, or references you can provide would be greatly appreciated.

Thank you in advance for your assistance!

I have already searched for relevant resources and documentation, but I couldn't find a clear example or tutorial specifically addressing this mixed copula model with the EM algorithm in R. Any guidance, code snippets, or references you can provide would be greatly appreciated.

Thank you in advance for your assistance!

Answer 1

Using the EM algorithm with a copula model is the same as the EM algorithm with any other models, such as a mixture of the Gaussian model.

If you do not prefer to do it manually, the scopula package in R include a mixture copula model that can fit only two types of copula functions at a time. Also, there is a MixCopula function in the copula package.

Also, you can use this code, where you can fit a mixture of Gaussian copulas manually:

  # Load required libraries
library(copula)
library(mvtnorm)

# Generate simulated data with a bivariate mixture distribution
set.seed(123)
n <- 1000
theta <- 0.3
mu1 <- c(0, 0)
mu2 <- c(3, 3)
sigma1 <- matrix(c(1, 0.5, 0.5, 2), ncol = 2)
sigma2 <- matrix(c(2, -0.5, -0.5, 1), ncol = 2)
u <- runif(n)
x <- ifelse(u < theta, MASS::mvrnorm(n, mu1, sigma1), MASS::mvrnorm(n, mu2, sigma2))

# Initialize the parameters
mu <- list(matrix(c(-1, -1), ncol = 2), matrix(c(1, 1), ncol = 2))
sigma <- list(matrix(c(1, 0.5, 0.5, 1), ncol = 2), matrix(c(2, -0.5, -0.5, 1), ncol = 2))
pi <- c(0.5, 0.5)

# Define the log-likelihood function
loglik <- function(x, mu, sigma, pi) {
  n <- dim(x)[1]
  likelihoods <- matrix(0, n, 2)
  for(i in 1:2) {
    likelihoods[,i] <- dmvnorm(x, mean = mu[[i]], sigma = sigma[[i]])
  }
  loglik <- sum(log(likelihoods %*% pi))
  return(loglik)
}

# Define the Gaussian copula function
gaussian_copula <- function(x, sigma) {
  n <- dim(x)[1]
  d <- dim(x)[2]
  rho <- cor(x)
  inv_sigma <- solve(sigma)
  u <- pnorm(x, mean = rep(0, d), sd = rep(1, d))
  z <- qnorm(u, mean = rep(0, d), sd = rep(1, d))
  phi <- dmvnorm(z, mean = rep(0, d), sigma = inv_sigma)
  copula <- pnorm(z, mean = rep(0, d), sd = rep(1, d))
  return(list(copula = copula, phi = phi))
}

# Run the EM algorithm
tolerance <- 1e-6
converged <- FALSE
while(!converged) {
  # E-step: compute the responsibilities
  likelihoods <- matrix(0, n, 2)
  copulas <- list(NULL, NULL)
  for(i in 1:2) {
    copulas[[i]] <- gaussian_copula(x, sigma[[i]])
    likelihoods[,i] <- copulas[[i]]$phi * pi[i]
  }
  responsibilities <- likelihoods / rowSums(likelihoods)
  
  # M-step: update the parameters
  N <- rowSums(responsibilities)
  pi <- N / n
  for(i in 1:2) {
    mu[[i]] <- colSums(responsibilities[,i] * x) / N[i]
    z <- qnorm(copulas[[i]]$copula, mean = rep(0, 2), sd = rep(1, 2))
    sigma[[i]] <- solve(var(z) * (1 - copulas[[i]]$copula)^2)
  }
  
  # Check for convergence
  old_loglik <- loglik(x, mu, sigma, pi)
  new_loglik <- loglik(x, mu, sigma, pi)
  if(abs(new_loglik - old_loglik) < tolerance) {
    converged <- TRUE
  }
}

# Output the estimated parameters
cat("mu1 =", mu[[1]], "mu2 =", mu[[2]], "sigma1 =", sigma[[1]], "sigma2 =", sigma[[2]], "pi1 =", pi[1], "pi2 =", pi[2])