R/calculate_lsq.R

In GRIDCOPULA: Bivariate Copula Functions Based on Regular Grid

calculate_lsq <- function(U, k, m) {
  Qm <- count.grid(U, k, m)
  Am <- (m*k) * (Qm / sum(Qm))
  
  # limSolve package
  Quad.m <- diag(k*m)
  Quad.v <- as.vector(Am)
  Eq.m <- matrix(0, nrow=k+m, ncol=k*m)
  for(i in 1:m) {
    Eq.m[i, ((i-1)*k+1):(i*k)] <- 1
  }
  for(i in 1:k) {
    Eq.m[(m+i), seq(i, m*k, by=k)] <- 1
  }
  Eq.v <- c(rep(k, m), rep(m, k))
  Eq.m <- Eq.m[-nrow(Eq.m), ]
  Eq.v <- Eq.v[-length(Eq.v)]
  Ineq.m <- rbind(diag(k*m), -diag(k*m))
  Ineq.v <- c(rep(0,k*m), -rep(min(k,m),k*m))
  value <- limSolve::lsei(A=Quad.m, B=Quad.v, E=Eq.m, F=Eq.v, G=Ineq.m, H=Ineq.v, 
                          type=1)
  Dm <- matrix(value$X, nrow=k, ncol=m)
  
  return(list(Density=Dm, Quantity=Qm, m=m, k=k))
}