RobustCovariance.Rmd

---
title: "R: Robust Covariance"
author: "Max Chen 2018"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Library

Load required packages

```{r lib}
# install.packages("mvtnorm")
# install.packages("robust")
# install.packages("car")
# install.packages("pracma")
# install.packages("ggplot2")
suppressMessages(library(mvtnorm))
suppressMessages(library(robust))
suppressMessages(library(car))
suppressMessages(library(pracma))
suppressMessages(library(ggplot2))
```

## Data

Generate sample data

```{r dat}
rndsample = rmvnorm(100,mean=c(3,5),sigma=matrix(c(1,0.5,0.5,2),2,2))
colnames(rndsample) = c("X1","X2")
df = data.frame(rndsample)
par(pty="s")
plot(df,asp=1,xlim=c(-2,8),ylim=c(0,10))
```

## Calculation

Calculate robust covariance estimator

```{r cal}
robfit = covRob(df)
print(robfit$center)
print(robfit$cov)
```

Transform back to unit distribution

```{r unit}
decov = svd(robfit$cov)
matS = diag(sqrt(decov$d))
matR = decov$u
matC = robfit$center
df_unit = as.matrix(sweep(df,2,matC,"-")) %*% solve(matR) %*% solve(matS)
par(pty="s")
plot(df_unit,asp=1)
```

## Visualize

Plot confidence ellipse

```{r plt}
nDim = 2
cumP = 0.95
radiProb = sqrt(qchisq(cumP,nDim))
dist = sqrt(rowSums(df_unit^2))
radiEvlp = max(dist)
par(pty="s")
plot(df,asp=1,xlim=c(-2,8),ylim=c(0,10))
ellipse(robfit$center,robfit$cov,radiProb)
ellipse(robfit$center,robfit$cov,radiEvlp,lty="dashed")
```

Plot density function

```{r pdf}
XY = meshgrid(seq(-2,8,length=100),seq(0,10,length=100))
XY = data.frame(as.vector(XY$X),as.vector(XY$Y))
colnames(XY) = c("X1","X2")
XY$pdf = dmvnorm(XY,mean=robfit$center,sigma=robfit$cov)
ggplot(XY,aes(X1,X2)) + 
    geom_raster(aes(fill=pdf),interpolate=TRUE) + 
    scale_fill_gradientn(colours=c("white","blue"))
```