Package {mvdpd}

Robust estimation for multivariate normal model using DPD based methods

Description

Robust estimation for multivariate normal model using multivariate (ordinary), composite and componentwise Density Power Divergence (DPD) miminization.

Usage

  mvnormDPD(x, beta = 0.3,
    method = c("multivariate", "composite", "componentwise"),
    initial = "MAD", ...)

Arguments

Details

Iterative algorithms are used to compute the estimates for all methods if beta >0. Additional algorithmic control parameters common to all methods are tol representing the tolerance limit so that the algorithm is said to converge if the total aboslute errors across all parameters are less than or equal to it (default value 1e-3), and max.it denoting the maximum number of iterations to be tried by the algorithm (default 1000).

For method="multivariate" and method="composite", additional control parameters are w.lower, w.upper. Weights smaller than w.lower are set to zero while weights larger than w.upper are set to one.

For method="composite" at the initial step a univariate filter can be applied setting univariate.filter=TRUE. In this case a univariate DPD with beta equals to univariate.beta (default 0.99) is computed for each column, data are then standardized and values larger than q are set to NA.

Value

An object of class mvdpd, that for this model is a list with components

Author(s)

A. Ghosh, C. Agostinelli, and A. Basu.

References

A. Basu, I.R. Harris, N.L. Hjort, and M.C. Jones. Robust and efficient estimation by minimizing a density power divergence. Biometrika, 85(3), 549-559, 1998.

A. Basu, A. Ghosh, and L. Pardo. Statistical Inference based on the Density Power Divergence: The Robustness Perspective. CRC Press, 2025.

S. Chakraborty, A. Basu, and A. Ghosh. A componentwise estimation procedure for multivariate location and scatter: Robustness, efficiency and scalability. Journal of Multivariate Analysis, 105546, 2025.

A. Ghosh, C. Agostinelli, and A. Basu. A Composite Divergence Approach to Robust Multivariate Estimation under Cellwise and Casewise Contamination. arXiv:2608.18914, http://arxiv.org/abs/2608.18914, 2026.

Examples

  mu <- rep(0, 3)
  Sigma <- diag(3) * 0.5 + 0.5
  set.seed(123)
  X <- MASS::mvrnorm(1000, mu, Sigma)
  X[1:5, 1]  <- X[1:5, 1] + 5
  X[6:10, 2] <- X[6:10, 2] - 10
  X[12, 1:2] <- c(-4,8)
  colnames(X) <- c("X1","X2","X3")
  M <- mvnormDPD(X, beta=0.3, method="multivariate")
  C <- mvnormDPD(X, beta=0.3, method="composite")
  W <- mvnormDPD(X, beta=0.3, method="componentwise")
  rbind(M$mu, C$mu, W$mu)

The Multivariate Ordered Gamma-Generated Distribution

Description

Density and random generation for the multivariate ordered gamma-generated distribution with scale parameters given by delta and baseline rate equal to lambda. The bivariate case coincides with the Mckay bivariate gamma distribution.

Usage

dmvogamma(x, delta, lambda, beta=1, log = FALSE)
rmvogamma(n, delta, lambda)

Arguments

Details

delta and lambda are positive parameters. In dmvogamma the result is the density power to beta.

Value

dmvogamma gives the density, rmvogamma generates random deviates.

Author(s)

A. Ghosh, C. Agostinelli, and A. Basu.

References

N. Balakrishnan and M.M. Ristic. Multivariate families of gamma-generated distributions with finite or infinite support above or below the diagonal. Journal of Multivariate Analysis, 143,194-207, 2016.

See Also

Distributions for other standard distributions.

Examples

  set.seed(123)
  n <- 100  
  delta <- seq(3, 0.5, by=-0.6)
  p <- length(delta)  
  lambda <- 0.5
  x <- rmvogamma(n, delta, lambda)
  dmvogamma(x, delta, lambda)

Robust estimation for multivariate ordered gamma model using DPD based methods

Description

Robust estimation for multivariate gamma generated model (Balakrishnan and Ristic, 2016) for ordered positive data using composite Density Power Divergence (DPD) miminization.

Usage

  mvogammaDPD(x, beta = 0.3,
    method = c("multivariate", "composite"),
    initial = list(delta=NULL, lambda=NULL))

Arguments

Value

An object of class mvdpd, that for this model is a list with components

Author(s)

A. Ghosh, C. Agostinelli, and A. Basu.

References

N. Balakrishnan and M.M. Ristic. Multivariate families of gamma-generated distributions with finite or infinite support above or below the diagonal. Journal of Multivariate Analysis, 143,194-207, 2016.

A. Basu, A. Ghosh, and L. Pardo. Statistical Inference based on the Density Power Divergence: The Robustness Perspective. CRC Press, 2025.

A. Ghosh, C. Agostinelli, and A. Basu. A Composite Divergence Approach to Robust Multivariate Estimation under Cellwise and Casewise Contamination. arXiv2608.18914, http://arxiv.org/abs/2608.18914, 2026.

See Also

mvogamma for the description of the multivariate ordered gamma generated distributions. See also twocenturies50 for an application of the method to a real dataset about Ultra Marathon Races

Examples

  set.seed(123)
  n <- 80
  delta <- seq(3, 0.5, by=-0.6)
  p <- length(delta)  
  lambda <- 0.5
  x <- rmvogamma(n=n, delta=delta, lambda=lambda)
  mvogammaDPD(x, initial=list(delta=delta, lambda=lambda), beta=0.3)

Maximum Likelihood estimation for multivariate ordered gamma model

Description

Maximum Likelihood for a multivariate gamma generated model (Balakrishnan and Ristic, 2016) for ordered positive data.

Usage

  mvogammaML(x, method=c("multivariate", "composite"),
    initial = list(delta=NULL, lambda=NULL))

Arguments

Value

An object of class mvdpd, that for this model is a list with components

Author(s)

A. Ghosh, C. Agostinelli, and A. Basu.

References

N. Balakrishnan and M.M. Ristic. Multivariate families of gamma-generated distributions with finite or infinite support above or below the diagonal. Journal of Multivariate Analysis, 143,194-207, 2016.

See Also

mvogamma for the description of the multivariate ordered gamma generated distributions and mvogammaDPD for a robust method of estimation

Examples

  set.seed(123)
  n <- 80
  delta <- seq(3, 0.5, by=-0.6)
  p <- length(delta)  
  lambda <- 0.5
  x <- rmvogamma(n=n, delta=delta, lambda=lambda)
  mvogammaML(x, initial=list(delta=delta, lambda=lambda))

Summary Method for "mvdpd" Objects

Description

Summary method for R object of class "mvdpd".

Usage

  ## S3 method for class 'mvdpd'
summary(object, se = c("sandwich", "exact"),
    null.hypothesis=c(rep(0,p), rep(1,p), rep(0, p*(p-1)/2)),
    print.outliers = FALSE, ...)

Arguments

Details

For pbjects with model=="normal" identification of cell, couple and row outliers are also reported by comparing the Mahalanobis distance (univariate, bivariate and row-wise) with the corresponding chi-quare distribution.

Value

summary(object) returns an object of S3 class "summary.mvdpd", basically a list that contains all the elements from the input object of class mvdpd and the following elements

Author(s)

A. Ghosh, C. Agostinelli, and A. Basu.

References

A. Basu, I.R. Harris, N.L. Hjort, and M.C. Jones. Robust and efficient estimation by minimizing a density power divergence. Biometrika, 85(3), 549-559, 1998.

A. Basu, A. Ghosh, and L. Pardo. Statistical Inference based on the Density Power Divergence: The Robustness Perspective. CRC Press, 2025.

S. Chakraborty, A. Basu, and A. Ghosh. A componentwise estimation procedure for multivariate location and scatter: Robustness, efficiency and scalability. Journal of Multivariate Analysis, 105546, 2025.

A. Ghosh, C. Agostinelli, and A. Basu. A Composite Divergence Approach to Robust Multivariate Estimation under Cellwise and Casewise Contamination. arXiv2608.18914, http://arxiv.org/abs/2608.18914, 2026.

Examples

  mu <- rep(0, 3)
  Sigma <- diag(3) * 0.5 + 0.5
  set.seed(123)
  X <- MASS::mvrnorm(1000, mu, Sigma)
  colnames(X) <- c("X1","X2","X3")
  M <- mvnormDPD(X, beta=0.3, method="multivariate")
  summary(M, se="sandwich")
  summary(M, se="exact", print.outliers=TRUE)

Two Centuries of Ultra Marathon Races (Male 50-miles races)

Description

This dataset was obtained from the publicly available data on the website https://www.kaggle.com/datasets/fatihyavuzz/two-centuries-of-um-races, which contains ace results from ultramarathon events held worldwide over approximately two centuries.

We restrict to mens' 50-mile races. For each race event (identified by the event name and date), finishing times were converted to hours and the ten fastest finishers were retained, producing one ordered $10$-dimensional observation per race. Events with fewer than ten recorded finishers or with tied finishing times among the top ten were discarded, yielding a collection of complete ordered observations of 4683 races.

Usage

data("twocenturies50")

Format

A numeric matrix with 4683 rows and 10 columns. Each row contains the best ten recorded finishers in a given race.

Source

https://www.kaggle.com/datasets/fatihyavuzz/two-centuries-of-um-races

Examples

  data(twocenturies50)
  ordered_data <- twocenturies50[order(twocenturies50[,1]),]
  n <- nrow(twocenturies50)
  p <- ncol(twocenturies50)

  # Best 200 races in term of the first finisher (clean data?)
  best0 <- ordered_data[1:200,]
  best0 <- best0/3600
  # Maximum likelihood
  ml <- mvogammaML(best0, initial=list(delta=c(18, rep(0.5, 9)),
    lambda=5))
  ml
  

  
  # Maximum composite likelihood (minimum CDPD with beta=0)
  cml <- mvogammaDPD(x=best0, beta=0, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  # Minimum composite DPD
  cdpd1 <- mvogammaDPD(x=best0, beta=0.1, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd2 <- mvogammaDPD(x=best0, beta=0.2, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd3 <- mvogammaDPD(x=best0, beta=0.3, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd5 <- mvogammaDPD(x=best0, beta=0.5, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  # Minimum (multivariate) DPD
  mdpd1 <- mvogammaDPD(x=best0, beta=0.1, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd2 <- mvogammaDPD(x=best0, beta=0.2, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd3 <- mvogammaDPD(x=best0, beta=0.3, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd5 <- mvogammaDPD(x=best0, beta=0.5, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))

  results0.multivariate <- cbind(
    c(ml$delta, ml$lambda),
    c(mdpd1$delta, mdpd1$lambda),
    c(mdpd2$delta, mdpd2$lambda),
    c(mdpd3$delta, mdpd3$lambda),
    c(mdpd5$delta, mdpd5$lambda)
  )
  
  results0.composite <- cbind(
    c(cml$delta, cml$lambda),
    c(cdpd1$delta, cdpd1$lambda),
    c(cdpd2$delta, cdpd2$lambda),
    c(cdpd3$delta, cdpd3$lambda),
    c(cdpd5$delta, cdpd5$lambda)
  )

  colnames(results0.multivariate) <- c("ML", "DPD(0.1)", "DPD(0.2)",
    "DPD(0.3)", "DPD(0.5)")
  colnames(results0.composite) <- c("CML", "CDPD(0.1)", "CDPD(0.2)",
    "CDPD(0.3)", "CDPD(0.5)")
  rownames(results0.multivariate) <- c(paste0("delta", 1:10), "lambda")
  rownames(results0.composite) <- rownames(results0.multivariate) 
  
  # Casewise contaminated data
  best1 <- ordered_data[c(1:180, (n-19):n),]
  best1 <- best1/3600
  # Maximum likelihood
  ml <- mvogammaML(best1, initial=list(delta=c(18, rep(0.5, 9)),
    lambda=5))
  # Maximum composite likelihood (minimum CDPD with beta=0)
  cml <- mvogammaDPD(x=best1, beta=0, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  # Minimum composite DPD
  cdpd1 <- mvogammaDPD(x=best1, beta=0.1, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd2 <- mvogammaDPD(x=best1, beta=0.2, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd3 <- mvogammaDPD(x=best1, beta=0.3, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd5 <- mvogammaDPD(x=best1, beta=0.5, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  # Minimum (multivariate) DPD
  mdpd1 <- mvogammaDPD(x=best1, beta=0.1, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd2 <- mvogammaDPD(x=best1, beta=0.2, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd3 <- mvogammaDPD(x=best1, beta=0.3, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd5 <- mvogammaDPD(x=best1, beta=0.5, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))

  results1.multivariate <- cbind(
    c(ml$delta, ml$lambda),
    c(mdpd1$delta, mdpd1$lambda),
    c(mdpd2$delta, mdpd2$lambda),
    c(mdpd3$delta, mdpd3$lambda),
    c(mdpd5$delta, mdpd5$lambda)
  )
  
  results1.composite <- cbind(
    c(cml$delta, cml$lambda),
    c(cdpd1$delta, cdpd1$lambda),
    c(cdpd2$delta, cdpd2$lambda),
    c(cdpd3$delta, cdpd3$lambda),
    c(cdpd5$delta, cdpd5$lambda)
  )

  colnames(results1.multivariate) <- c("ML", "DPD(0.1)", "DPD(0.2)",
    "DPD(0.3)", "DPD(0.5)")
  colnames(results1.composite) <- c("CML", "CDPD(0.1)", "CDPD(0.2)",
    "CDPD(0.3)", "CDPD(0.5)")
  rownames(results1.multivariate) <- c(paste0("delta", 1:10), "lambda")
  rownames(results1.composite) <- rownames(results1.multivariate) 

  # Cellwise contaminated data
  set.seed(1234)
  best2 <- best0
  for (i in 1:200) {
    nn <- rbinom(n=1, size=10, prob=0.1)
    if (nn){
      pos <- sample(1:p, size=nn)
      best2[i,pos] <- runif(1, 0, 3)*best0[i,pos]
      best2[i,] <- sort(best2[i,]) 
    }
  }
  best2 <- best2/3600
  # Maximum likelihood
  ml <- mvogammaML(best2, initial=list(delta=c(18, rep(0.5, 9)),
    lambda=5))
  # Maximum composite likelihood (minimum CDPD with beta=0)
  cml <- mvogammaDPD(x=best2, beta=0, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  # Minimum composite DPD
  cdpd1 <- mvogammaDPD(x=best2, beta=0.1, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd2 <- mvogammaDPD(x=best2, beta=0.2, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd3 <- mvogammaDPD(x=best2, beta=0.3, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  cdpd5 <- mvogammaDPD(x=best2, beta=0.5, method="composite",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  # Minimum (multivariate) DPD
  mdpd1 <- mvogammaDPD(x=best2, beta=0.1, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd2 <- mvogammaDPD(x=best2, beta=0.2, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd3 <- mvogammaDPD(x=best2, beta=0.3, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))
  mdpd5 <- mvogammaDPD(x=best2, beta=0.5, method="multivariate",
    initial=list(delta=ml$delta, lambda=ml$lambda))

  results2.multivariate <- cbind(
    c(ml$delta, ml$lambda),
    c(mdpd1$delta, mdpd1$lambda),
    c(mdpd2$delta, mdpd2$lambda),
    c(mdpd3$delta, mdpd3$lambda),
    c(mdpd5$delta, mdpd5$lambda)
  )
  
  results2.composite <- cbind(
    c(cml$delta, cml$lambda),
    c(cdpd1$delta, cdpd1$lambda),
    c(cdpd2$delta, cdpd2$lambda),
    c(cdpd3$delta, cdpd3$lambda),
    c(cdpd5$delta, cdpd5$lambda)
  )

  colnames(results2.multivariate) <- c("ML", "DPD(0.1)", "DPD(0.2)",
    "DPD(0.3)", "DPD(0.5)")
  colnames(results2.composite) <- c("CML", "CDPD(0.1)", "CDPD(0.2)",
    "CDPD(0.3)", "CDPD(0.5)")
  rownames(results2.multivariate) <- c(paste0("delta", 1:10), "lambda")
  rownames(results2.composite) <- rownames(results2.multivariate)