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()`).

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()`).

plotAlcohol$eigen_df_gg

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

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()`).

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()`).

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()`).

plotMilk$eigen_df_gg

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

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()`).

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()`).

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()`).

plotBushfire$eigen_df_gg

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

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()`).

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