library(expm)
## Loading required package: Matrix
##
## Attaching package: 'expm'
## The following object is masked from 'package:Matrix':
##
## expm
A <- matrix(c(25,-2,4,-2,4,1,4,1,9), 3, 3, byrow=TRUE)
A
## [,1] [,2] [,3]
## [1,] 25 -2 4
## [2,] -2 4 1
## [3,] 4 1 9
E <- eigen(A)
E$values
## [1] 26.078452 8.495796 3.425752
Inv <- solve(A)
Inv
## [,1] [,2] [,3]
## [1,] 0.04611331 0.02898551 -0.02371542
## [2,] 0.02898551 0.27536232 -0.04347826
## [3,] -0.02371542 -0.04347826 0.12648221
B <- sqrtm(Inv)
B%*%B
## [,1] [,2] [,3]
## [1,] 0.04611331 0.02898551 -0.02371542
## [2,] 0.02898551 0.27536232 -0.04347826
## [3,] -0.02371542 -0.04347826 0.12648221
lm <- diag(E$values)
P <- E$vectors
P%*%lm%*%t(P)
## [,1] [,2] [,3]
## [1,] 25 -2 4
## [2,] -2 4 1
## [3,] 4 1 9
SQ <- sqrtm(A)
SQ%*%SQ
## [,1] [,2] [,3]
## [1,] 25 -2 4
## [2,] -2 4 1
## [3,] 4 1 9
A <- matrix(c(9,-2,-2,6), 2, 2, byrow=TRUE)
A
## [,1] [,2]
## [1,] 9 -2
## [2,] -2 6
E <- eigen(A)
E$values
## [1] 10 5
E$vectors
## [,1] [,2]
## [1,] -0.8944272 -0.4472136
## [2,] 0.4472136 -0.8944272
E$vectors[,1]%*%t(E$vectors[,1])
## [,1] [,2]
## [1,] 0.8 -0.4
## [2,] -0.4 0.2
E$vectors[,2]%*%t(E$vectors[,2])
## [,1] [,2]
## [1,] 0.2 0.4
## [2,] 0.4 0.8
E$vectors%*%diag(E$values) %*%t(E$vectors)
## [,1] [,2]
## [1,] 9 -2
## [2,] -2 6
\[ 10\begin{bmatrix} 0.8 & -0.4 \\ -0.4 & 0.2 \end{bmatrix} + 5\begin{bmatrix} 0.2 & 0.4 \\ 0.4 & 0.8 \end{bmatrix} \]
A <- matrix(c(25,-2,4,-2,4,1,4,1,9), 3, 3, byrow=TRUE)
diag(A)
## [1] 25 4 9
sqrtm(diag(diag(A)))
## [,1] [,2] [,3]
## [1,] 5 0 0
## [2,] 0 2 0
## [3,] 0 0 3
cov2cor(A)
## [,1] [,2] [,3]
## [1,] 1.0000000 -0.2000000 0.2666667
## [2,] -0.2000000 1.0000000 0.1666667
## [3,] 0.2666667 0.1666667 1.0000000
covr <- diag(c(3,3,3,3))
covr
## [,1] [,2] [,3] [,4]
## [1,] 3 0 0 0
## [2,] 0 3 0 0
## [3,] 0 0 3 0
## [4,] 0 0 0 3
A <- matrix(c(1,-1,0,0,1,1,-2,0,1,1,1,-3), 3, 4, byrow=TRUE)
A
## [,1] [,2] [,3] [,4]
## [1,] 1 -1 0 0
## [2,] 1 1 -2 0
## [3,] 1 1 1 -3
mu <- matrix(c(3,2,-2,0), 4, 1, byrow=TRUE)
mu
## [,1]
## [1,] 3
## [2,] 2
## [3,] -2
## [4,] 0
A%*%mu
## [,1]
## [1,] 1
## [2,] 9
## [3,] 3
A%*%covr%*%t(A)
## [,1] [,2] [,3]
## [1,] 6 0 0
## [2,] 0 18 0
## [3,] 0 0 36
Each of the linear combinations have zero covariances, as seen by the Cov(Ax) we found which only has numbers on the diagonal. Those being the variances of each Ax.
\[ \Sigma^{-1/2}\mathbf{X} \sim \mathcal{N}(\mathbf{0}, \Sigma^{-1/2}\Sigma\Sigma^{-1/2}) = \mathcal{N}(\mathbf{0}, \mathbf{Q}\Lambda^{-1/2}\mathbf{Q}^\top\mathbf{Q}\Lambda\mathbf{Q}^\top\mathbf{Q}\Lambda^{-1/2}\mathbf{Q}^\top) \] \[ \qquad\qquad\qquad\qquad\qquad\qquad\qquad = \mathcal{N}(\mathbf{0}, \mathbf{Q}\Lambda^{-1/2}\Lambda\Lambda^{-1/2}\mathbf{Q}^\top) = \mathcal{N}(\mathbf{0}, \mathbf{Q}\mathbf{Q}^\top) = \mathcal{N}(\mathbf{0}, \mathbf{I}_p) \] Therefore, \[ X^\top\Sigma^{-1}X = (\Sigma^{-1/2}X)^\top(\Sigma^{-1/2}X) \sim \chi^2_p \]
E[X] = E[3\(X_1\)-2\(X_2\)+\(X_3\)] = 3E[\(X_1\)]-2E[\(X_2\)]+E[\(X_3\)] = (3(2)) - (2(-3)) + 1 = 6 + 6 + 1 = 13
Cov(X) = Cov(3\(X_1\)-2\(X_2\)+\(X_3\)) = 9Cov(\(X_1\)) + 4Cov(\(X_2\)) + Cov(\(X_3\)) - 6Cov(\(X_1\),\(X_2\)) - 2Cov(\(X_2\),\(X_3\)) = 9(1) + 4(3) - 6(1) - 2(2) + 2 = 5
So, X follows a normal distribution with mean = 13 and variance = 5.
library(MASS)
Sig <- matrix(c(1,1,1,1,3,2,1,2,2), 3, 3, byrow=TRUE)
Sig
## [,1] [,2] [,3]
## [1,] 1 1 1
## [2,] 1 3 2
## [3,] 1 2 2
mun <- matrix(c(2,-3,1), 1, 3, byrow=TRUE)
mun
## [,1] [,2] [,3]
## [1,] 2 -3 1
mvrnorm(n=100, mu=mun, Sigma=Sig)
## [,1] [,2] [,3]
## [1,] 2.3769072 -3.8022287 0.67927117
## [2,] 2.8121345 -2.3858826 1.51808193
## [3,] 0.5427505 -3.7667466 -0.74777526
## [4,] 0.9022711 -5.7630559 -0.36282428
## [5,] 4.6683065 0.9457601 4.29174438
## [6,] -0.4873659 -5.0468999 -1.55689576
## [7,] 2.4096913 -2.7348166 1.49763075
## [8,] 4.1634869 -0.1602428 3.72943398
## [9,] 3.3449334 -1.4738125 2.85329741
## [10,] 0.9267578 -1.3854292 1.32781759
## [11,] 2.4441280 -4.1141280 1.47153080
## [12,] -0.6034607 -6.2829383 -0.58108332
## [13,] 2.4977472 -2.1784651 2.62538691
## [14,] 2.6269220 -2.2263561 1.86415685
## [15,] 3.9574999 -1.6842188 1.94392455
## [16,] 0.6234990 -3.8359612 0.23888569
## [17,] 3.1792919 -4.4348679 0.86160346
## [18,] 2.8606617 -4.1429598 0.41745803
## [19,] 3.5367355 -0.9397996 2.85773491
## [20,] 1.0645574 -3.6158340 0.33188889
## [21,] 2.1175881 -1.9400208 0.59095521
## [22,] 1.1734328 -5.0658330 -1.38885026
## [23,] 0.5652790 -5.0311559 -0.42608539
## [24,] 2.5976596 -1.9807694 1.46141722
## [25,] 1.9486238 -2.4598526 0.85651282
## [26,] 0.8760408 -4.5872967 -1.06832266
## [27,] 2.5238365 -1.7994215 0.76059917
## [28,] 4.7224709 -0.3780642 4.45923887
## [29,] 3.8017517 -3.2560666 0.76080742
## [30,] 0.1646429 -3.6154488 0.67072960
## [31,] 2.1291733 -4.4085951 0.02760084
## [32,] 3.3692766 -2.5945001 1.07490447
## [33,] 1.5591462 -4.3762859 0.49703596
## [34,] 1.0498793 -2.4309113 0.33956877
## [35,] 2.5042414 -4.0087583 2.20717500
## [36,] 2.9411905 -3.0490442 0.66785853
## [37,] 0.5471877 -3.8440516 -0.10052473
## [38,] 3.0239777 -4.1686970 0.43682398
## [39,] 2.3926222 -3.1504494 1.02567819
## [40,] 0.1755333 -6.4047311 -2.31980020
## [41,] 3.7873176 -1.9813476 2.35051659
## [42,] 2.1324222 -0.9420295 3.03829055
## [43,] 2.5083254 -3.7050782 2.08872463
## [44,] 2.3965690 -1.2026981 2.65750893
## [45,] 2.5366661 -2.5469773 1.39173783
## [46,] 3.4986215 -0.7297347 3.58607632
## [47,] 2.4070826 -3.7882981 2.22505777
## [48,] 1.9138290 -3.7368599 1.08384508
## [49,] 3.8720630 -4.3127000 2.02858036
## [50,] 1.2093045 -3.7356670 0.79006550
## [51,] 1.6207370 -5.4621482 -0.74344449
## [52,] 2.5113580 -3.1875717 2.79153366
## [53,] 2.8321805 0.1098860 2.54417673
## [54,] 1.1324378 -3.3145793 0.62202849
## [55,] 3.1259417 -1.4050760 3.18759422
## [56,] 1.7275307 -2.4493007 1.53818187
## [57,] 1.0891497 -5.5106388 -1.03305975
## [58,] 2.4247322 -2.3293467 0.36653226
## [59,] 2.1604353 -2.5252917 1.45493972
## [60,] 1.3882703 -3.0494442 1.26369914
## [61,] 1.8043284 -5.6757220 -0.62831410
## [62,] 1.4476699 -3.0358717 0.72899185
## [63,] 0.8268250 -2.6752629 0.25040627
## [64,] 3.2276639 -2.4404583 2.25035247
## [65,] 0.7268224 -2.2877365 -0.46879968
## [66,] 2.5236475 -4.3799829 1.05576466
## [67,] 4.4768919 2.0085720 6.01299092
## [68,] 2.8230409 -0.7773724 3.10499315
## [69,] 1.8190918 -5.9821129 -0.15826919
## [70,] 0.8699412 -4.7755101 0.14356103
## [71,] 2.4068765 -2.9359900 1.94174969
## [72,] 3.0532191 -0.2437066 3.54908590
## [73,] 2.8929926 -2.7365869 2.07077706
## [74,] 2.7010246 -3.7379791 1.50371553
## [75,] 3.0316001 -3.2115182 1.60382039
## [76,] -0.2483000 -3.5343431 -1.48646756
## [77,] 2.8031077 -2.0647857 1.42393060
## [78,] 2.0419343 -3.1931870 0.22103038
## [79,] 1.3479309 -3.9552413 0.57202598
## [80,] 1.5912353 -3.8988121 1.06224014
## [81,] 1.4817359 -4.0180357 1.34200951
## [82,] 0.3262046 -5.1625191 -0.18204045
## [83,] 0.9996796 -5.6607590 -0.73617238
## [84,] 1.1199066 -2.7652038 0.21597840
## [85,] 1.2747338 -3.3015216 0.21150347
## [86,] 0.6468893 -2.5329120 0.49612811
## [87,] 1.4730682 -4.7612501 -0.20217078
## [88,] 1.6208512 -3.5294566 0.84644599
## [89,] 2.4457240 -1.8760914 2.17164569
## [90,] -0.1801747 -7.0107378 -2.80008337
## [91,] 2.0757575 -2.6484666 2.04887370
## [92,] 1.1096680 -5.5526943 -1.41228705
## [93,] 1.7127258 -3.7796973 1.43912545
## [94,] 1.0839349 -4.3600006 -0.01285655
## [95,] 2.3805904 -4.6713308 0.97473439
## [96,] 0.7790854 -5.4298371 -1.05264877
## [97,] 2.5307414 -1.7725410 1.73050726
## [98,] 2.0870693 -2.9676472 2.20968310
## [99,] 2.0486155 -1.9458749 2.30495447
## [100,] 1.6524113 -3.3668725 0.55022689
Sigm <- matrix(c(1,-2,0,-2,5,0,0,0,2), 3, 3, byrow=TRUE)
Sigm
## [,1] [,2] [,3]
## [1,] 1 -2 0
## [2,] -2 5 0
## [3,] 0 0 2
muy <- matrix(c(-3,1,4), 3, 1, byrow=TRUE)
muy
## [,1]
## [1,] -3
## [2,] 1
## [3,] 4
#Only ones with covariance are X(1) and X(2), which is -2
Cov(\(\frac{X_1+X_2}{2}\), \(X_3\)) = \(\frac{1}{2}\)Cov(\(X_1\),\(X_3\))+\(\frac{1}{2}\)Cov(\(X_2\),\(X_3\)) = \(\frac{1}{2}\)(0) + \(\frac{1}{2}\)(0) = 0 Since their covariance is 0, they are independent.
Cov(\(X_2\), \(X_2\) - \(X_1\) - \(X_3\)) = Var(\(X_2\)) - Cov(\(X_1\),\(X_2\)) - Cov(\(X_2\),\(X_3\)) = 5 - (-2) - 0 = 7 Since their covariance is 7, they are not independent.
x5 <- c(12,9,5,8,8,12,12,21,11,13,10,12,18,11,8,9,7,16,13,9, 14,7,13,5,10,7,11,7,9,7,10,12,8,10,6,9,6,13,9,8,11,6)
x6 <- c(8,5,6,15,10,12,15,14,11,9,3,7,10,7,10,10,7,4,2,5,4,6,11,2,23,6,11,10,8,2,7,8,4,24,9,10,12,18,25,6,14,5)
comb <- cbind(x5,x6)
Xf <- as.matrix(comb)
Xf
## x5 x6
## [1,] 12 8
## [2,] 9 5
## [3,] 5 6
## [4,] 8 15
## [5,] 8 10
## [6,] 12 12
## [7,] 12 15
## [8,] 21 14
## [9,] 11 11
## [10,] 13 9
## [11,] 10 3
## [12,] 12 7
## [13,] 18 10
## [14,] 11 7
## [15,] 8 10
## [16,] 9 10
## [17,] 7 7
## [18,] 16 4
## [19,] 13 2
## [20,] 9 5
## [21,] 14 4
## [22,] 7 6
## [23,] 13 11
## [24,] 5 2
## [25,] 10 23
## [26,] 7 6
## [27,] 11 11
## [28,] 7 10
## [29,] 9 8
## [30,] 7 2
## [31,] 10 7
## [32,] 12 8
## [33,] 8 4
## [34,] 10 24
## [35,] 6 9
## [36,] 9 10
## [37,] 6 12
## [38,] 13 18
## [39,] 9 25
## [40,] 8 6
## [41,] 11 14
## [42,] 6 5
muf <- colMeans(Xf)
muf
## x5 x6
## 10.047619 9.404762
sigmaf <- cov(Xf)
sigmaf
## x5 x6
## x5 11.363531 3.126597
## x6 3.126597 30.978513
#Distances
d <- NULL
for (i in 1:42){
d <- c(d, t(Xf[i,]-muf) %*%solve(sigmaf) %*%(Xf[i,]-muf))
}
d
## [1] 0.4606524 0.6592206 2.3770610 1.6282902 0.4135364 0.4760726
## [7] 1.1848895 10.6391792 0.1388339 0.8162468 1.3566301 0.6228096
## [13] 5.6494392 0.3159498 0.4135364 0.1224973 0.8987982 4.7646873
## [19] 3.0089122 0.6592206 2.7741416 1.0360061 0.7874152 3.4437748
## [25] 6.1488606 1.0360061 0.1388339 0.8856041 0.1379719 2.2488867
## [31] 0.1901188 0.4606524 1.1471939 7.0857237 1.4584229 0.1224973
## [37] 1.8984708 2.7782596 8.4730649 0.6370218 0.7032485 1.8013611
#Proportion
n <- nrow(Xf)
p <- ((1:n)-0.5)/n
p
## [1] 0.01190476 0.03571429 0.05952381 0.08333333 0.10714286 0.13095238
## [7] 0.15476190 0.17857143 0.20238095 0.22619048 0.25000000 0.27380952
## [13] 0.29761905 0.32142857 0.34523810 0.36904762 0.39285714 0.41666667
## [19] 0.44047619 0.46428571 0.48809524 0.51190476 0.53571429 0.55952381
## [25] 0.58333333 0.60714286 0.63095238 0.65476190 0.67857143 0.70238095
## [31] 0.72619048 0.75000000 0.77380952 0.79761905 0.82142857 0.84523810
## [37] 0.86904762 0.89285714 0.91666667 0.94047619 0.96428571 0.98809524
#Chi-Square Plot
chsq <- qchisq(p,df=2, lower.tail=TRUE)
plot(sort(d),chsq)
X <- as.matrix(iris[,1:4])
X
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## [1,] 5.1 3.5 1.4 0.2
## [2,] 4.9 3.0 1.4 0.2
## [3,] 4.7 3.2 1.3 0.2
## [4,] 4.6 3.1 1.5 0.2
## [5,] 5.0 3.6 1.4 0.2
## [6,] 5.4 3.9 1.7 0.4
## [7,] 4.6 3.4 1.4 0.3
## [8,] 5.0 3.4 1.5 0.2
## [9,] 4.4 2.9 1.4 0.2
## [10,] 4.9 3.1 1.5 0.1
## [11,] 5.4 3.7 1.5 0.2
## [12,] 4.8 3.4 1.6 0.2
## [13,] 4.8 3.0 1.4 0.1
## [14,] 4.3 3.0 1.1 0.1
## [15,] 5.8 4.0 1.2 0.2
## [16,] 5.7 4.4 1.5 0.4
## [17,] 5.4 3.9 1.3 0.4
## [18,] 5.1 3.5 1.4 0.3
## [19,] 5.7 3.8 1.7 0.3
## [20,] 5.1 3.8 1.5 0.3
## [21,] 5.4 3.4 1.7 0.2
## [22,] 5.1 3.7 1.5 0.4
## [23,] 4.6 3.6 1.0 0.2
## [24,] 5.1 3.3 1.7 0.5
## [25,] 4.8 3.4 1.9 0.2
## [26,] 5.0 3.0 1.6 0.2
## [27,] 5.0 3.4 1.6 0.4
## [28,] 5.2 3.5 1.5 0.2
## [29,] 5.2 3.4 1.4 0.2
## [30,] 4.7 3.2 1.6 0.2
## [31,] 4.8 3.1 1.6 0.2
## [32,] 5.4 3.4 1.5 0.4
## [33,] 5.2 4.1 1.5 0.1
## [34,] 5.5 4.2 1.4 0.2
## [35,] 4.9 3.1 1.5 0.2
## [36,] 5.0 3.2 1.2 0.2
## [37,] 5.5 3.5 1.3 0.2
## [38,] 4.9 3.6 1.4 0.1
## [39,] 4.4 3.0 1.3 0.2
## [40,] 5.1 3.4 1.5 0.2
## [41,] 5.0 3.5 1.3 0.3
## [42,] 4.5 2.3 1.3 0.3
## [43,] 4.4 3.2 1.3 0.2
## [44,] 5.0 3.5 1.6 0.6
## [45,] 5.1 3.8 1.9 0.4
## [46,] 4.8 3.0 1.4 0.3
## [47,] 5.1 3.8 1.6 0.2
## [48,] 4.6 3.2 1.4 0.2
## [49,] 5.3 3.7 1.5 0.2
## [50,] 5.0 3.3 1.4 0.2
## [51,] 7.0 3.2 4.7 1.4
## [52,] 6.4 3.2 4.5 1.5
## [53,] 6.9 3.1 4.9 1.5
## [54,] 5.5 2.3 4.0 1.3
## [55,] 6.5 2.8 4.6 1.5
## [56,] 5.7 2.8 4.5 1.3
## [57,] 6.3 3.3 4.7 1.6
## [58,] 4.9 2.4 3.3 1.0
## [59,] 6.6 2.9 4.6 1.3
## [60,] 5.2 2.7 3.9 1.4
## [61,] 5.0 2.0 3.5 1.0
## [62,] 5.9 3.0 4.2 1.5
## [63,] 6.0 2.2 4.0 1.0
## [64,] 6.1 2.9 4.7 1.4
## [65,] 5.6 2.9 3.6 1.3
## [66,] 6.7 3.1 4.4 1.4
## [67,] 5.6 3.0 4.5 1.5
## [68,] 5.8 2.7 4.1 1.0
## [69,] 6.2 2.2 4.5 1.5
## [70,] 5.6 2.5 3.9 1.1
## [71,] 5.9 3.2 4.8 1.8
## [72,] 6.1 2.8 4.0 1.3
## [73,] 6.3 2.5 4.9 1.5
## [74,] 6.1 2.8 4.7 1.2
## [75,] 6.4 2.9 4.3 1.3
## [76,] 6.6 3.0 4.4 1.4
## [77,] 6.8 2.8 4.8 1.4
## [78,] 6.7 3.0 5.0 1.7
## [79,] 6.0 2.9 4.5 1.5
## [80,] 5.7 2.6 3.5 1.0
## [81,] 5.5 2.4 3.8 1.1
## [82,] 5.5 2.4 3.7 1.0
## [83,] 5.8 2.7 3.9 1.2
## [84,] 6.0 2.7 5.1 1.6
## [85,] 5.4 3.0 4.5 1.5
## [86,] 6.0 3.4 4.5 1.6
## [87,] 6.7 3.1 4.7 1.5
## [88,] 6.3 2.3 4.4 1.3
## [89,] 5.6 3.0 4.1 1.3
## [90,] 5.5 2.5 4.0 1.3
## [91,] 5.5 2.6 4.4 1.2
## [92,] 6.1 3.0 4.6 1.4
## [93,] 5.8 2.6 4.0 1.2
## [94,] 5.0 2.3 3.3 1.0
## [95,] 5.6 2.7 4.2 1.3
## [96,] 5.7 3.0 4.2 1.2
## [97,] 5.7 2.9 4.2 1.3
## [98,] 6.2 2.9 4.3 1.3
## [99,] 5.1 2.5 3.0 1.1
## [100,] 5.7 2.8 4.1 1.3
## [101,] 6.3 3.3 6.0 2.5
## [102,] 5.8 2.7 5.1 1.9
## [103,] 7.1 3.0 5.9 2.1
## [104,] 6.3 2.9 5.6 1.8
## [105,] 6.5 3.0 5.8 2.2
## [106,] 7.6 3.0 6.6 2.1
## [107,] 4.9 2.5 4.5 1.7
## [108,] 7.3 2.9 6.3 1.8
## [109,] 6.7 2.5 5.8 1.8
## [110,] 7.2 3.6 6.1 2.5
## [111,] 6.5 3.2 5.1 2.0
## [112,] 6.4 2.7 5.3 1.9
## [113,] 6.8 3.0 5.5 2.1
## [114,] 5.7 2.5 5.0 2.0
## [115,] 5.8 2.8 5.1 2.4
## [116,] 6.4 3.2 5.3 2.3
## [117,] 6.5 3.0 5.5 1.8
## [118,] 7.7 3.8 6.7 2.2
## [119,] 7.7 2.6 6.9 2.3
## [120,] 6.0 2.2 5.0 1.5
## [121,] 6.9 3.2 5.7 2.3
## [122,] 5.6 2.8 4.9 2.0
## [123,] 7.7 2.8 6.7 2.0
## [124,] 6.3 2.7 4.9 1.8
## [125,] 6.7 3.3 5.7 2.1
## [126,] 7.2 3.2 6.0 1.8
## [127,] 6.2 2.8 4.8 1.8
## [128,] 6.1 3.0 4.9 1.8
## [129,] 6.4 2.8 5.6 2.1
## [130,] 7.2 3.0 5.8 1.6
## [131,] 7.4 2.8 6.1 1.9
## [132,] 7.9 3.8 6.4 2.0
## [133,] 6.4 2.8 5.6 2.2
## [134,] 6.3 2.8 5.1 1.5
## [135,] 6.1 2.6 5.6 1.4
## [136,] 7.7 3.0 6.1 2.3
## [137,] 6.3 3.4 5.6 2.4
## [138,] 6.4 3.1 5.5 1.8
## [139,] 6.0 3.0 4.8 1.8
## [140,] 6.9 3.1 5.4 2.1
## [141,] 6.7 3.1 5.6 2.4
## [142,] 6.9 3.1 5.1 2.3
## [143,] 5.8 2.7 5.1 1.9
## [144,] 6.8 3.2 5.9 2.3
## [145,] 6.7 3.3 5.7 2.5
## [146,] 6.7 3.0 5.2 2.3
## [147,] 6.3 2.5 5.0 1.9
## [148,] 6.5 3.0 5.2 2.0
## [149,] 6.2 3.4 5.4 2.3
## [150,] 5.9 3.0 5.1 1.8
mu <- colMeans(X)
mu
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 5.843333 3.057333 3.758000 1.199333
sigma <- cov(X)
sigma
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 0.6856935 -0.0424340 1.2743154 0.5162707
## Sepal.Width -0.0424340 0.1899794 -0.3296564 -0.1216394
## Petal.Length 1.2743154 -0.3296564 3.1162779 1.2956094
## Petal.Width 0.5162707 -0.1216394 1.2956094 0.5810063
d <- NULL
for (i in 1:4){
d <- c(d, t(X[i,]-mu) %*%solve(sigma) %*%(X[i,]-mu))
}
d
## [1] 2.134468 2.849119 2.081339 2.452382
p <- ((1:4)-0.5)/4
p
## [1] 0.125 0.375 0.625 0.875
chsq <- qchisq(p,df=4, lower.tail=TRUE)
plot(sort(d),chsq)
abline(a=0, b=1, lwd=2)
#Outlier Detection
cut_point <- qchisq(0.05,df=4,lower.tail=FALSE)
d[(d>cut_point)]
## numeric(0)
The first 4 appear to roughly follow a multivariate normal distribution, although the plot is a bit misleading as it is zoomed too far in. There are no outliers as seen by the numeric(0).
#Centering
XX <- sweep(X, 2, mu, FUN="-")
SVD <- svd(XX)
SVD$d
## [1] 25.099960 6.013147 3.413681 1.884524
X2 <- XX%*%SVD$v
zapsmall(cov(X2))
## [,1] [,2] [,3] [,4]
## [1,] 4.228242 0.000000 0.00000 0.000000
## [2,] 0.000000 0.242671 0.00000 0.000000
## [3,] 0.000000 0.000000 0.07821 0.000000
## [4,] 0.000000 0.000000 0.00000 0.023835
We can observe that because there is only numbers on the diagonal, that there are no covariances and therefore the linear combinations are independent.
E <- eigen(sigma)
diag(E$values)
## [,1] [,2] [,3] [,4]
## [1,] 4.228242 0.0000000 0.0000000 0.00000000
## [2,] 0.000000 0.2426707 0.0000000 0.00000000
## [3,] 0.000000 0.0000000 0.0782095 0.00000000
## [4,] 0.000000 0.0000000 0.0000000 0.02383509
The eigenvalues of the covariance matrix of X are therefore equal to linear combinations variances.