partial.mahalanobis: Partial Square Mahalanobis Distance

In GSE: Robust Estimation in the Presence of Cellwise and Casewise Contamination and Missing Data

Partial Square Mahalanobis Distance

Description

Computes the partial square Mahalanobis distance for all observations in x. Let \mathbf{x} = (x_{i1},...,x_{ip})' be a p-dimensional random vector and \mathbf{u} = (u_{i1},...,u_{ip})' be a p-dimensional vectors of zeros and ones indicating which entry is missing: 0 as missing and 1 as observed. Then partial mahalanobis distance is given by:

d(\mathbf{x, u, m, Σ}) = (\mathbf{x^{(u)}} - \mathbf{m^{(u)}})' (\mathbf{Σ^{(u)}})^{-1}(\mathbf{x^{(u)}} - \mathbf{m^{(u)}}).

With no missing data, this function is equivalent to mahalanobis distance.

Usage

partial.mahalanobis(x, mu, Sigma)

Arguments

x	a matrix or data frame. May contain missing values, but cannot contain columns with completely missing entries.
mu	location estimate
Sigma	scatter estimate. Must be positive definite

Value

An S4 object of class CovRobMiss-class. The output S4 object contains the following slots:

mu	Estimated location. Can be accessed via getLocation.
S	Estimated scatter matrix. Can be accessed via getScatter.
pmd	Squared partial Mahalanobis distances. Can be accessed via getDist.
pmd.adj	Adjusted squared partial Mahalanobis distances. Can be accessed via getDistAdj.
pu	Dimension of the observed entries for each case. Can be accessed via getDim.
call	Object of class "language". Not meant to be accessed.
x	Input data matrix. Not meant to be accessed.
p	Column dimension of input data matrix. Not meant to be accessed.
estimator	Character string of the name of the estimator used. Not meant to be accessed.

Author(s)

Andy Leung [e-mail]

Examples

## Not run: 
## suppose we would like to compute pmd for an MLE
x <- matrix(rnorm(1000),100,10)
U <- matrix(rbinom(1000,1,0.1),100,10)
x <- x * ifelse(U==1,NA,1)
## compute MLE (i.e. EM in this case)
res <- CovEM(x)
## compute pmd
res.pmd <- partial.mahalanobis(x, mu=getLocation(res), S=getScatter(res))
summary(res.pmd)
plot(res.pmd, which="index")

## End(Not run)

GSE documentation built on Dec. 28, 2022, 1:31 a.m.