confidence ellipses

t.r

> library(MASS)    # mvrnorm
> library(nortest) # ad.test
> library(car)     # dataEllipse
> library(emdbook) # HPDregionplot
> library(MVN)     # mardiaTest,hzTest,roystonTest
> ukff_tmp <- read.csv('sampled_ukff.csv')
> ukff <- ukff_tmp[complete.cases(ukff_tmp),]

> colMeans(ukff)
         X0        X0.1 
-0.02565391 -1.21415781 
> var(ukff)
            X0      X0.1
X0    11.18782 -10.70098
X0.1 -10.70098  22.45431
> summary(ukff)
       X0                 X0.1         
 Min.   :-12.11505   Min.   :-17.8708  
 1st Qu.: -2.49540   1st Qu.: -3.9409  
 Median :  0.00000   Median : -0.9968  
 Mean   : -0.02565   Mean   : -1.2142  
 3rd Qu.:  2.71686   3rd Qu.:  1.9694  
 Max.   :  8.44793   Max.   : 13.7895  
> # Simulate process for multivariate normal distribution with the same mean and covariance matrix
> ran <- mvrnorm(n=1000,colMeans(ukff), var(ukff))
> colMeans(ran)
        X0       X0.1 
-0.1007181 -1.3161636 
> var(ran)
            X0      X0.1
X0    11.20242 -11.44782
X0.1 -11.44782  23.40849
> summary(ran)
       X0                X0.1        
 Min.   :-12.1232   Min.   :-15.357  
 1st Qu.: -2.2786   1st Qu.: -4.556  
 Median : -0.1277   Median : -1.625  
 Mean   : -0.1007   Mean   : -1.316  
 3rd Qu.:  2.1433   3rd Qu.:  1.792  
 Max.   :  9.1622   Max.   : 13.627  

> dataEllipse(ran[,1], ran[,2], levels=c(.9973,.9545,.6827), main='test')
> plot(ran)
> HPDregionplot(ran, prob = c(0.9973, 0.9545, 0.6827), col=c("salmon", "lightblue", "lightgreen"), lwd=3, add=TRUE)
> shapiro.test(ran[,1])

    Shapiro-Wilk normality test

data:  ran[, 1]
W = 0.99694, p-value = 0.05202

> shapiro.test(ran[,2])

    Shapiro-Wilk normality test

data:  ran[, 2]
W = 0.99697, p-value = 0.05434
> ad.test(ran[,1])

    Anderson-Darling normality test

data:  ran[, 1]
A = 0.69744, p-value = 0.06843

> ad.test(ran[,2])

    Anderson-Darling normality test

data:  ran[, 2]
A = 0.86168, p-value = 0.0269

> qqnorm(ran[,1])
> qqline(ran[,1])
> qqnorm(ran[,2])
> qqline(ran[,2])
> mardiaTest(ran, qqplot=TRUE)
   Mardia's Multivariate Normality Test 
--------------------------------------- 
   data : ran 

   g1p            : 0.03080181 
   chi.skew       : 5.133636 
   p.value.skew   : 0.2738584 

   g2p            : 7.848846 
   z.kurtosis     : -0.5974903 
   p.value.kurt   : 0.5501801 

   chi.small.skew : 5.159345 
   p.value.small  : 0.271335 

   Result         : Data are multivariate normal. 
--------------------------------------- 

> roystonTest(ran)
  Royston's Multivariate Normality Test 
--------------------------------------------- 
  data : ran 

  H       : 5.202394 
  p-value : 0.0425277 

  Result  : Data are not multivariate normal. 
--------------------------------------------- 

> hzTest(ran)
  Henze-Zirkler's Multivariate Normality Test 
--------------------------------------------- 
  data : ran 

  HZ      : 0.8546706 
  p-value : 0.2935467 

  Result  : Data are multivariate normal. 
--------------------------------------------- 
# Ukff data
> plot(ukff)
> HPDregionplot(ukff, prob = c(0.9973, 0.9545, 0.6827), col=c("salmon", "lightblue", "lightgreen"), lwd=3, add=TRUE)
> dev.new()
> dataEllipse(ukff[,1], ukff[,2], levels=c(.9973,.9545,.6827), main='test')
> shapiro.test(ukff[,1])

    Shapiro-Wilk normality test

data:  ukff[, 1]
W = 0.98933, p-value = 1.827e-06

> shapiro.test(ukff[,2])

    Shapiro-Wilk normality test

data:  ukff[, 2] 
W = 0.9945, p-value = 0.001394

> ad.test(ukff[,1])

    Anderson-Darling normality test

data:  ukff[, 1]
A = 3.425, p-value = 1.44e-08

> ad.test(ukff[,2])

    Anderson-Darling normality test

data:  ukff[, 2]
A = 1.734, p-value = 0.0001921

> qqnorm(ukff[,1])
> qqline(ukff[,1])
> qqnorm(ukff[,2])
> qqline(ukff[,2])
> dev.new()
> mardiaTest(ukff, qqplot=TRUE)
   Mardia's Multivariate Normality Test
--------------------------------------- 
   data : ukff 

   g1p            : 1.062354 
   chi.skew       : 170.5079 
   p.value.skew   : 8.136973e-36 

   g2p            : 9.889232 
   z.kurtosis     : 7.328386 
   p.value.kurt   : 2.329248e-13 

   chi.small.skew : 171.3946 
   p.value.small  : 5.249651e-36 

   Result          : Data are not multivariate normal.
--------------------------------------- 

> hzTest(ukff)
  Henze-Zirkler's Multivariate Normality Test
--------------------------------------------- 
  data : ukff

  HZ      : 7.131201
  p-value : 0

  Result  : Data are not multivariate normal.
--------------------------------------------- 

> roystonTest(ukff)
  Royston's Multivariate Normality Test 
--------------------------------------------- 
  data : ukff 

  H       : 25.74584
  p-value : 1.210226e-06

  Result  : Data are not multivariate normal.
---------------------------------------------