-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy path01-09-Stat-Methods-Eigen-Distance.Rmd
69 lines (53 loc) · 1.95 KB
/
01-09-Stat-Methods-Eigen-Distance.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
## Eigen Values and Statistical Distance {-}
```{r a25, warning=FALSE, message=FALSE, comment=NA}
## Eigen values and vectors are used to describe a positie definate covariance matrix.
(Cov = matrix(c(1, 0, 0, 2), nrow = 2))
eig = eigen(Cov)
## eigen values
(lambda = eig$values)
## eigen vectors
(ee = eig$vectors)
## Spectural decompsotion allows you to reconstruct a matrix using only the eigen
## values and vectors
lambda[1] * ee[,1] %*% t(ee[,1]) + lambda[2] * ee[,2] %*% t(ee[,2])
## Straight line distance (Euclidean) vs Statistical Distance
## Straight line distance to the origin using point(1, 1)
sqrt((1 - 0)^2 + (1 - 0)^2)
# Statistical Distance to the origin
sqrt((1 - 0)^2/1 + (1 - 0)^2/5)
## Generating Multivariate Normal Data
library(mvtnorm)
## set up parameters, 2 means use the covariance matrix from earlier
mu = c(10, 20)
## generate a large dataset
set.seed(1000)
X = rmvnorm(1000, mu, Cov)
## Correlation of X, should be close to 0
cor(X)
## Calculate the distance between each point and the means
distance = c()
for (i in 1:nrow(X)) {
x = t(X[i, ] - colMeans(X)) %*% solve(cov(X)) %*% (X[i, ] - colMeans(X))
distance = c(distance, x)
}
## distances
head(distance)
## What is the distance that captures 95% of all points generated from the distribution?
(critical.value = qchisq(.95, 2))
## What is the proportion of points that fall within this distance?
## As n increases the proportion should converge on 5%
length(which((distance - critical.value) > 0))/length(distance)
## plot a 95% confidence ellipse for the generated data
library(plotrix)
plot(x = c(7,13), y = c(15,24), type = "n",
xlab = expression(x[1]), ylab = expression(x[2]),
main = "95% Confidence Ellipsoid")
points(X[, 1], X[, 2])
abline(h = 20, v = 10, lty = 2, lwd = 2, col = "red")
draw.ellipse(10, 20,
sqrt(critical.value * lambda[1]),
sqrt(critical.value * lambda[2]),
## convert radians to degrees
angle = acos(abs(ee[1,1])) * 57.2957795,
border = 1, lwd = 2)
```