mvnormDPD examples

A. Ghosh, C. Agostinelli and A. Basu

2026-08-24

Introduction

This file contains some examples of the use of the function mvnormDPD. It reproduces all the figures and tables in section C “Performances under benchmark Gaussian datasets” of A. Ghosh, C. Agostinelli and A. Basu (2026) A Composite Divergence Approach to Robust Multivariate Estimation under Cellwise and Casewise Contamination, arXiv:2608.18914, https://arxiv.org/abs/2608.18914.

library("mvdpd")
library("cellWise")
library("robustbase")
library("dplyr")

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library("ggplot2")
library("knitr")

Functions and wrappers

loo.analysis <- function(X, method) {
  p <- ncol(X)
  var_names <- colnames(X)
  eigen_res <- list()
  mu_res <- list()
  md_res <- list()
  ## Leave-One-Out
  for (i in 1:nrow(X)) {
    Xb <- X[-i,]
    res <- method(Xb)
    mu_res[[i]] <- data.frame(
      LeaveOut = i,
      Method = res$method,
      Parameter = var_names,
      Estimate = res$mu
    )
    rownames(mu_res[[i]]) <- var_names
  
    eigen_res[[i]] <- data.frame(
      LeaveOut = i,
      Method = res$method,
      Parameter = paste0("e", 1:p),
      Eigenvalues = eigen(res$Sigma)$values
    )
    rownames(eigen_res[[i]]) <- paste0("e", 1:p)
  
    md_res[[i]] <- data.frame(
      LeaveOut = i,
      Method = res$method,
      Observation = paste0("Obs", seq_len(nrow(X))),
      Distance = sqrt(mahalanobis(X,
        center = res$mu,
        cov = res$Sigma))
    )
  }  
  res <- list(mu=mu_res, eigen=eigen_res, mahalanobis=md_res)
  return(res)
}

ML <- function(X, ...) {
  list(mu=colMeans(X), Sigma=cov(X), method="ML/MCL")
}

MCD <- function(X, ...) { 
  mcd <- covMcd(X, ...)
  list(mu=mcd$center, Sigma=mcd$cov, method="MCD")
}

CELLMCD <- function(X, ...) { 
  cellmcd <- cellMCD(X, checkPars=list(silent=TRUE), ...)
  list(mu=cellmcd$mu, Sigma=cellmcd$S, method="CellMCD")
}

MDPD <- function(X, beta, ...) { 
  mdpd <- mvnormDPD(X, beta=beta, method="multivariate", ...)
  list(mu=mdpd$mu, Sigma=mdpd$Sigma, method=paste0("MDPD(", beta, ")"))
}

CDPD <- function(X, beta, ...) { 
  cdpd <- mvnormDPD(X, beta=beta, method="composite", ...)
  list(mu=cdpd$mu, Sigma=cdpd$Sigma, method=paste0("CDPD(", beta, ")"))
}

perform.analysis <- function(X) {  
  resML <- loo.analysis(X, method=ML)
  resMCD <- loo.analysis(X, method=MCD)
  resCELLMCD <- loo.analysis(X, method=CELLMCD)
  resMDPD1 <- loo.analysis(X, method=function(X) MDPD(X, beta=0.1))
  resMDPD3 <- loo.analysis(X, method=function(X) MDPD(X, beta=0.3))
  resMDPD5 <- loo.analysis(X, method=function(X) MDPD(X, beta=0.5))
  resCDPD1 <- loo.analysis(X, method=function(X) CDPD(X, beta=0.1))
  resCDPD3 <- loo.analysis(X, method=function(X) CDPD(X, beta=0.3))
  resCDPD5 <- loo.analysis(X, method=function(X) CDPD(X, beta=0.5))
  mu_df <- bind_rows(resML$mu, resMCD$mu, resCELLMCD$mu, 
    resMDPD1$mu, resMDPD3$mu, resMDPD5$mu,
    resCDPD1$mu, resCDPD3$mu, resCDPD5$mu)
  eigen_df <- bind_rows(resML$eigen, resMCD$eigen, resCELLMCD$eigen, 
    resMDPD1$eigen, resMDPD3$eigen, resMDPD5$eigen,
    resCDPD1$eigen, resCDPD3$eigen, resCDPD5$eigen)

  md_df <- bind_rows(resML$mahalanobis, 
    resMCD$mahalanobis, resCELLMCD$mahalanobis, 
    resMDPD1$mahalanobis, resMDPD3$mahalanobis, resMDPD5$mahalanobis,
    resCDPD1$mahalanobis, resCDPD3$mahalanobis, resCDPD5$mahalanobis)
  
  mu_var <- mu_df %>%
    group_by(Method, Parameter) %>%
    summarise(
      Variance = var(Estimate),
      .groups = "drop"
    )

  eigen_var <- eigen_df %>%
    group_by(Method, Parameter) %>%
    summarise(
      Variance = var(Eigenvalues),
      .groups = "drop"
    )

  md_var <- md_df %>%
    group_by(Method, Observation) %>%
    summarise(
      Variance = var(Distance),
      .groups = "drop"
    )
  res <- list(mu_df=mu_df, mu_var=mu_var,
    eigen_df=eigen_df, eigen_var=eigen_var,
    md_df=md_df,  md_var=md_var)
  return(res)
}

plot.results <- function(object) {
  # Boxplot of the location estimates
  mu_df_gg <- ggplot(object$mu_df, aes(x=Method, y=Estimate, fill=Parameter)) +
    geom_boxplot(position = position_dodge(0.8), width = 0.7) +
    scale_x_discrete(name="Method", limits=c("CellMCD", "ML/MCL", "CDPD(0.1)", "CDPD(0.3)", "CDPD(0.5)")) +
    scale_y_continuous(name=expression(paste("Leave-one-out ",hat(mu)))) + 
    theme(
      axis.text = element_text(size = 12, face="bold"),
      axis.title.x = element_text(size = 14, face="bold"),
      axis.title.y = element_text(size = 16, face="bold"),
      legend.text = element_text(size = 12, face="bold"),
      legend.title = element_text(size = 14, face="bold")
    )

  # Boxplot of variances of the location estimates
  mu_var_gg <- ggplot(object$mu_var, aes(x = Method, y = Variance)) +
    geom_boxplot(fill = "grey85", width = 0.6) +
    geom_point(aes(color = Parameter), size = 3, 
      position = position_jitter(width = 0.08)) +
    scale_x_discrete(name="Method",
      limits=c("CellMCD", "ML/MCL", "CDPD(0.1)", "CDPD(0.3)", "CDPD(0.5)")) +
    scale_y_continuous(name=expression(paste("Variances of leave-one-out ",
      hat(mu)))) + 
    theme(
      axis.text = element_text(size = 12, face="bold"),
      axis.title.x = element_text(size = 14, face="bold"),
      axis.title.y = element_text(size = 16, face="bold"),
      legend.text = element_text(size = 12, face="bold"),
      legend.title = element_text(size = 14, face="bold")
    )

  # Boxplot of the eigenvalues of Scatter estimates
  eigen_df_gg <- ggplot(object$eigen_df, aes(x=Method,
    y=Eigenvalues, fill=Parameter)) +
    geom_boxplot(position = position_dodge(0.8), width = 0.7) +
    scale_x_discrete(name="Method", limits=c("CellMCD", "ML/MCL",
      "CDPD(0.1)", "CDPD(0.3)", "CDPD(0.5)")) +
    scale_y_continuous(name=expression(paste("Leave-one-out eigenvalues of ",
      hat(Sigma)))) + 
    theme(
      axis.text = element_text(size = 12, face="bold"),
      axis.title.x = element_text(size = 14, face="bold"),
      axis.title.y = element_text(size = 16, face="bold"),
      legend.text = element_text(size = 12, face="bold"),
      legend.title = element_text(size = 14, face="bold")
    )

  # Boxplot of variances of the eigenvalues of Scatter estimates
  eigen_var_gg <- ggplot(object$eigen_var, aes(x = Method, y = Variance)) +
    geom_boxplot(fill = "grey85", width = 0.6) +
    geom_point(aes(color = Parameter), size = 3, 
      position = position_jitter(width = 0.08)) +
    scale_x_discrete(name="Method", limits=c("CellMCD",
      "ML/MCL", "CDPD(0.1)", "CDPD(0.3)", "CDPD(0.5)")) +
    scale_y_continuous(name=expression(paste("Variances of leave-one-out 
      eigenvalues of ",hat(Sigma)))) + 
    theme(
      axis.text = element_text(size = 12, face="bold"),
      axis.title.x = element_text(size = 14, face="bold"),
      axis.title.y = element_text(size = 16, face="bold"),
      legend.text = element_text(size = 12, face="bold"),
      legend.title = element_text(size = 14, face="bold")
    )
  res <- list(mu_df_gg=mu_df_gg, eigen_df_gg=eigen_df_gg,
    mu_var_gg=mu_var_gg, eigen_var_gg=eigen_var_gg)
  return(res)
}

Alcohol (n=44, d=7)

data(alcohol)
X <- as.matrix(alcohol)
X <- transfo(X)$Y

##  
##  The input data has 44 rows and 7 columns.

resAlcohol <- perform.analysis(X)
plotAlcohol <- plot.results(resAlcohol)

Figure S1

plotAlcohol$mu_df_gg

## Warning: Removed 1232 rows containing missing values or values outside the scale range
## (`stat_boxplot()`).
plot of chunk unnamed-chunk-12

plot of chunk unnamed-chunk-12

plotAlcohol$mu_var_gg

## Warning: Removed 28 rows containing missing values or values outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 28 rows containing missing values or values outside the scale range
## (`geom_point()`).
plot of chunk unnamed-chunk-13

plot of chunk unnamed-chunk-13

plotAlcohol$eigen_df_gg

## Warning: Removed 1232 rows containing missing values or values outside the scale range
## (`stat_boxplot()`).
plot of chunk unnamed-chunk-14

plot of chunk unnamed-chunk-14

plotAlcohol$eigen_var_gg

## Warning: Removed 28 rows containing missing values or values outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 28 rows containing missing values or values outside the scale range
## (`geom_point()`).
plot of chunk unnamed-chunk-15

plot of chunk unnamed-chunk-15

Milk (n=68, d=8)

data(milk)
X <- as.matrix(milk)
X <- transfo(X)$Y

##  
##  The input data has 86 rows and 8 columns.

resMilk <- perform.analysis(X)
plotMilk <- plot.results(resMilk)

Figure S2

plotMilk$mu_df_gg

## Warning: Removed 2752 rows containing missing values or values outside the scale range
## (`stat_boxplot()`).
plot of chunk unnamed-chunk-17

plot of chunk unnamed-chunk-17

plotMilk$mu_var_gg

## Warning: Removed 32 rows containing missing values or values outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 32 rows containing missing values or values outside the scale range
## (`geom_point()`).
plot of chunk unnamed-chunk-18

plot of chunk unnamed-chunk-18

plotMilk$eigen_df_gg

## Warning: Removed 2752 rows containing missing values or values outside the scale range
## (`stat_boxplot()`).
plot of chunk unnamed-chunk-19

plot of chunk unnamed-chunk-19

plotMilk$eigen_var_gg

## Warning: Removed 32 rows containing missing values or values outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 32 rows containing missing values or values outside the scale range
## (`geom_point()`).
plot of chunk unnamed-chunk-20

plot of chunk unnamed-chunk-20

Bushfire (n=38, d=5)

data(bushfire)
X <- as.matrix(bushfire)
X <- transfo(X)$Y

##  
##  The input data has 38 rows and 5 columns.

resBushfire <- perform.analysis(X)
plotBushfire <- plot.results(resBushfire)

Figure S3

plotBushfire$mu_df_gg

## Warning: Removed 760 rows containing missing values or values outside the scale range
## (`stat_boxplot()`).
plot of chunk unnamed-chunk-22

plot of chunk unnamed-chunk-22

plotBushfire$mu_var_gg

## Warning: Removed 20 rows containing missing values or values outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 20 rows containing missing values or values outside the scale range
## (`geom_point()`).
plot of chunk unnamed-chunk-23

plot of chunk unnamed-chunk-23

plotBushfire$eigen_df_gg

## Warning: Removed 760 rows containing missing values or values outside the scale range
## (`stat_boxplot()`).
plot of chunk unnamed-chunk-24

plot of chunk unnamed-chunk-24

plotBushfire$eigen_var_gg

## Warning: Removed 20 rows containing missing values or values outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 20 rows containing missing values or values outside the scale range
## (`geom_point()`).
plot of chunk unnamed-chunk-25

plot of chunk unnamed-chunk-25

Toxicity (n=38, d=10)

Transformed data

data(toxicity)
X <- as.matrix(toxicity)
X <- transfo(X)$Y

##  
##  The input data has 38 rows and 10 columns.

betas <- c(0, 0.1, 0.3, 0.5)
d <- ncol(X)
pairs <- which(upper.tri(matrix(0, d, d)), arr.ind=TRUE)
result <-matrix(NA, nrow = d*(d+3)/2, ncol=length(betas))
for (i in seq_along(betas)){
  res <- mvnormDPD(X, betas[i], method = "composite")
  cor_matrix <- res$rho
  result[,i] <- c(res$mu, res$sigma,
    cor_matrix[upper.tri(cor_matrix, diag = FALSE)])
}
colnames(result) <- c("ML/MCL", "CDPD(0.1)", "CDPD(0.3)", "CDPD(0.5)")
rownames(result) <- c(paste0("mu_", 1:10), paste0("sigma2_", 1:10), paste0("rho_", pairs[,1], pairs[,2]))

Table S1

kable(result)
ML/MCL CDPD(0.1) CDPD(0.3) CDPD(0.5)
mu_1 -0.0779019 -0.0776080 -0.0852228 -0.0675296
mu_2 0.0000000 0.0164793 0.0580872 0.0952169
mu_3 0.0000000 0.0194761 0.0651009 0.1234987
mu_4 0.0000000 0.0243500 0.0771907 0.1321183
mu_5 -0.0705141 -0.0601355 -0.0394666 -0.0242888
mu_6 0.0000000 -0.0030182 -0.0076627 -0.0135830
mu_7 0.0000000 0.0082279 0.0223343 0.0205170
mu_8 0.0000000 0.0061019 0.0122583 0.0104286
mu_9 0.0000000 -0.0081548 -0.0255715 -0.0484908
mu_10 0.0000000 0.0185714 0.0608116 0.1080782
sigma2_1 0.8868514 1.1918660 1.2001601 1.1658771
sigma2_2 0.7368421 1.0066978 1.0232597 1.0092182
sigma2_3 0.7368421 1.0381067 1.1452420 1.2164303
sigma2_4 0.7368421 1.0210980 1.0971169 1.1507239
sigma2_5 0.8561501 1.1303067 1.0858797 1.0122745
sigma2_6 0.7368421 1.0259148 1.1060055 1.1706685
sigma2_7 0.7368421 0.9993916 1.0177457 1.0182493
sigma2_8 0.7368421 1.0254819 1.1116006 1.1860707
sigma2_9 0.7368421 1.0194204 1.0885973 1.1465985
sigma2_10 0.7368421 1.0571674 1.2131411 1.3600295
rho_12 0.8322388 0.8628886 0.9152557 0.9250386
rho_13 0.1131656 0.1196088 0.1370754 0.1562779
rho_23 0.2808114 0.2835286 0.2756927 0.2609005
rho_14 0.2952595 0.3218801 0.3858360 0.4427930
rho_24 0.4860845 0.4980747 0.5072281 0.5148794
rho_34 0.7236445 0.7414810 0.7722202 0.7998671
rho_15 0.4804663 0.4974810 0.5271693 0.5386596
rho_25 0.7615118 0.7569288 0.7388229 0.7193244
rho_35 0.6083289 0.6015518 0.5778396 0.5411152
rho_45 0.5926599 0.5786725 0.5443224 0.5066433
rho_16 0.7319192 0.7317449 0.7285395 0.7260053
rho_26 0.8571748 0.8566828 0.8368186 0.8218142
rho_36 0.0833179 0.0769015 0.0595573 0.0375573
rho_46 0.2272637 0.2210286 0.2121922 0.2111919
rho_56 0.6249486 0.6311714 0.6456317 0.6690204
rho_17 0.5647034 0.5885093 0.6302105 0.6790042
rho_27 0.8278655 0.8341521 0.8337057 0.8556621
rho_37 0.1745410 0.1651261 0.1398045 0.1112415
rho_47 0.2450245 0.2383051 0.2193476 0.2019986
rho_57 0.8196564 0.8235258 0.8184323 0.8091686
rho_67 0.8333421 0.8502673 0.8622851 0.8645118
rho_18 -0.1301424 -0.1226917 -0.1033393 -0.0792777
rho_28 -0.0359887 -0.0276313 -0.0104282 0.0251174
rho_38 -0.5998177 -0.6165466 -0.6389156 -0.6561084
rho_48 -0.5783522 -0.6171167 -0.6806671 -0.7294232
rho_58 -0.0083235 -0.0111258 -0.0168985 -0.0251884
rho_68 0.1943763 0.2116480 0.2427001 0.2652904
rho_78 0.4185087 0.4295602 0.4313967 0.4263750
rho_19 -0.3173434 -0.3047845 -0.2793955 -0.2451218
rho_29 -0.2886486 -0.2808820 -0.2347194 -0.1693658
rho_39 -0.6278474 -0.6317284 -0.6386160 -0.6525587
rho_49 -0.6453400 -0.6821225 -0.7577854 -0.8116444
rho_59 -0.2901517 -0.2722216 -0.2269667 -0.1763799
rho_69 0.0583660 0.0769415 0.1208629 0.1873380
rho_79 0.1406515 0.1651036 0.2103375 0.2433964
rho_89 0.8850675 0.8954339 0.8989220 0.9048678
rho_110 0.2277795 0.2233137 0.2190017 0.2146313
rho_210 0.3587042 0.3606548 0.3444764 0.3231322
rho_310 0.8612577 0.8696825 0.8780543 0.8789530
rho_410 0.8013970 0.8091754 0.8243439 0.8419843
rho_510 0.4966583 0.4971298 0.4958759 0.4985985
rho_610 0.0677998 0.0593673 0.0455917 0.0313605
rho_710 0.1012489 0.1016181 0.0967266 0.0908839
rho_810 -0.6822559 -0.7077478 -0.7347179 -0.7490808
rho_910 -0.6834668 -0.6871258 -0.6932972 -0.6997174

Original data (without transformation)

data(toxicity)
X <- as.matrix(toxicity)
betas <- c(0, 0.1, 0.3, 0.5)
d <- ncol(X)
pairs <- which(upper.tri(matrix(0, d, d)), arr.ind=TRUE)
result <-matrix(NA, nrow = d*(d+3)/2, ncol=length(betas))
for (i in seq_along(betas)){
  res <- mvnormDPD(X, betas[i], method = "composite")
  cor_matrix <- res$rho
  result[,i] <- c(res$mu, res$sigma,
    cor_matrix[upper.tri(cor_matrix, diag = FALSE)])
}
colnames(result) <- c("ML/MCL", "CDPD(0.1)", "CDPD(0.3)", "CDPD(0.5)")
rownames(result) <- c(paste0("mu_", 1:10), paste0("sigma2_", 1:10), paste0("rho_", pairs[,1], pairs[,2]))

Table S2

kable(result)
ML/MCL CDPD(0.1) CDPD(0.3) CDPD(0.5)
mu_1 -0.1557895 -0.1755423 -0.2067308 -0.2356793
mu_2 1.6668421 1.6348539 1.6644823 1.6129076
mu_3 0.6489342 0.6444936 0.9038378 0.8953938
mu_4 4.3442105 4.3855271 4.6532383 4.6340566
mu_5 17.1908263 17.1778033 17.2250379 17.2222028
mu_6 3.1191263 2.9800565 2.8412039 2.6743840
mu_7 34.2747368 33.8515598 33.4611542 33.1061037
mu_8 1.4457632 1.4490532 1.4488896 1.4530081
mu_9 38.1447368 38.2452858 32.9937452 33.2497057
mu_10 6.8075526 1.4589594 1.4566000 1.4568976
sigma2_1 0.1228241 0.1648587 0.1612789 0.1671877
sigma2_2 1.3394439 1.8511166 1.9776046 2.1979217
sigma2_3 0.1480713 0.2084506 0.0110871 0.0129568
sigma2_4 0.5302261 0.6078094 0.0633034 0.0703075
sigma2_5 0.3379609 0.4466209 0.4212268 0.4359295
sigma2_6 6.9743331 9.2346262 9.3477467 9.5926329
sigma2_7 114.7460143 151.4431450 148.8149667 155.4520783
sigma2_8 0.0003998 0.0005507 0.0006251 0.0006260
sigma2_9 69.6179164 87.8339734 2.0839588 2.5158479
sigma2_10 31.9733432 0.0002870 0.0004049 0.0004736
rho_12 0.8229998 0.8420820 0.9193282 0.9428273
rho_13 -0.0961598 -0.0990495 0.5117906 0.5249982
rho_23 0.1509435 0.1610507 0.4244994 0.4204528
rho_14 -0.0024717 0.0575391 0.5727432 0.6138614
rho_24 0.3692350 0.3608647 0.4810635 0.4884879
rho_34 0.5684439 0.5881605 0.8623857 0.8576459
rho_15 0.4348867 0.4529446 0.5219870 0.6006378
rho_25 0.7514938 0.7481595 0.7289631 0.7301583
rho_35 0.5595866 0.5669460 0.3834806 0.3956621
rho_45 0.6442904 0.6160204 0.3245595 0.3067424
rho_16 0.7225736 0.7364638 0.7433900 0.7535805
rho_26 0.8691874 0.8596619 0.8493307 0.8322227
rho_36 0.1144352 0.0905973 0.0555203 0.0248107
rho_46 0.1407676 0.1265196 0.0915442 0.0856556
rho_56 0.5896459 0.5852197 0.5896458 0.6041467
rho_17 0.5743438 0.5997138 0.6396515 0.6895283
rho_27 0.8261782 0.8285386 0.8310871 0.8338219
rho_37 0.2103176 0.1986373 0.0006108 -0.0022758
rho_47 0.2523277 0.2197944 -0.0341496 -0.0119319
rho_57 0.7600946 0.7652160 0.7677910 0.7817935
rho_67 0.8920588 0.8922827 0.8856304 0.8894553
rho_18 -0.0432016 -0.0738186 -0.1325072 -0.2293875
rho_28 -0.0044649 -0.0008877 0.0143843 -0.0479339
rho_38 -0.3695229 -0.3694220 -0.6553771 -0.6197531
rho_48 -0.3416130 -0.4258811 -0.7628937 -0.7947009
rho_58 -0.0051047 -0.0071040 -0.0199242 -0.0553162
rho_68 0.2452322 0.2129706 0.2087165 0.1227525
rho_78 0.4332253 0.4028166 0.3843388 0.3094479
rho_19 -0.4622954 -0.4468673 0.1305310 0.0057667
rho_29 -0.5045037 -0.4794080 0.1150527 0.0354298
rho_39 -0.2017107 -0.2118809 -0.7010589 -0.7171404
rho_49 -0.4978286 -0.5733830 -0.4688995 -0.5425924
rho_59 -0.3915145 -0.3413993 -0.2678894 -0.2898946
rho_69 -0.2061051 -0.1632387 0.3772019 0.2951604
rho_79 -0.0600534 -0.0015708 0.2462240 0.2019507
rho_89 0.7446842 0.7449957 0.6435060 0.6369507
rho_110 0.3979104 -0.2413620 -0.2300900 -0.2111631
rho_210 0.5067174 0.0279085 0.0040479 0.0203725
rho_310 0.6019403 0.1496060 -0.2810203 -0.2372763
rho_410 0.4916877 -0.0891624 -0.6077092 -0.5921160
rho_510 0.5575355 0.3309931 0.2882333 0.3018304
rho_610 0.2760859 0.1245816 0.0745546 0.0597530
rho_710 0.1906814 0.4462666 0.4026010 0.3748937
rho_810 -0.5628562 0.8275556 0.7752512 0.8035794
rho_910 -0.5080758 0.6227186 0.0289097 0.0766466