amv27 - Distribution for Sum of Random Vectors. Conditional Distribution is Multivariate Normal.

Help this channel to remain great! Donating to Patreon or Paypal can do this!

[xinglinli]

OMG, you channel is really a treasure for me.

[jeanselee]

in equation of mu^star, element 9 in sigma_xx covariance matrix is 3 ?

[statisticsmatt]

begin video 1

x<-seq(-6, 6, .1)
plot(x, dnorm(x, 0, .5), type='l', ylab="Density Value", main="Univariate Normal Densities")
lines(x, dnorm(x, mean=0, sd=1), col='red')
lines(x, dnorm(x, mean=2, sd=1), col='green')
lines(x, dnorm(x, mean=-4, sd=.7), col='blue')
legend("topright", c("Mean=0, StdDev=0.5", "Mean=0, StdDev=1.0", "Mean=2, StdDev=1.0",
"Mean=-4, StdDev=0.7"), lty=1, col=c('black', 'red', 'green', 'blue'))

h <- function(x, y, mx, my, sx, sy, p) { (1/(2* pi * sx*sy*sqrt((1-p^2))))*
- 2*p*((x-mx)/sx)*((y-my)/sy) +
(y-my)^2/sy^2})}

par(mfrow=c(2, 2))
mx<-0;my<-0;sx<-2;sy<-2;p<- 0.1;stp<-.5
S <- matrix(c(sx^2, p*sx*sy, p*sx*sy, sy^2), ncol=2)
x <- seq(-3*sx, 3*sx, stp)
y <- seq(-3*sy, 3*sy, stp)
z <- outer(x, y, function(x, y) h(x, y, mx, my, sx, sy, p))
contour(x, y, z, xlab="y1", ylab="y2", main="Contour Plot of Bivariate Normal")
legend("bottomright", paste("Generalized Var = |Cov(y)| =", det(S)), bty='n')
persp(x, y, z, xlab="y1", ylab="y2", theta = 45, phi = 30, main="Density Plot. Correlation=0.1")

p<- 0.8
S <- matrix(c(sx^2, p*sx*sy, p*sx*sy, sy^2), ncol=2)
x <- seq(-3*sx, 3*sx, stp)
y <- seq(-3*sy, 3*sy, stp)
z <- outer(x, y, function(x, y) h(x, y, mx, my, sx, sy, p))
contour(x, y, z, xlab="y1", ylab="y2", main="Contour Plot of Bivariate Normal")
legend("bottomright", paste("Generalized Var = |Cov(y)| =", det(S)), bty='n')
persp(x, y, z, xlab="y1", ylab="y2", theta = -20, phi = 30, main="Density Plot. Correlation=0.8")
par(mfrow=c(1, 1))

end video 1

begin video 2

library(MASS) # load library MASS
n <- 1000
mu <- c(10, 5, 7, 9, 20)
Sig <- matrix(c(
5, 4, 3, 2, 1,
4, 5, 4, 3, 2,
3, 4, 5, 4, 3,
2, 3, 4, 5, 4,
1, 2, 3, 4, 5), ncol = 5)
Y <- mvrnorm(n, mu, Sig) # sample mvn data
tail(Y) ## last 6 observations

a <- c(1, 0, -1, 0, 0) # comparison between 1st and 3rd components
mu1 <- as.numeric(t(a)%*%mu) # mu1=3
Sig1 <- as.numeric(t(a)%*%Sig%*%a) # Sig1=4
z <- as.vector(t(a)%*%t(Y))
hist(z, probability = TRUE, xlab="a'Y",
main="Multivariate Normal to Univariate Normal")
lines(seq(-4, 12, .01), dnorm(seq(-4, 12, .01), mu1, sqrt(Sig1)), col='red', lwd=2)
legend("topleft", c("Mean=3", "Var=4"), title="Theoretical")
legend("topright", c(paste0("Mean=", round(mean(z), 1)),
paste0("Var=", round(var(z), 1))), title="Sample Stats")

library(MASS) # load library MASS
n <- 1000
mu <- c(10, 5, 7, 9, 20)
Sig <- matrix(c(
5, 4, 3, 2, 1,
4, 5, 4, 3, 2,
3, 4, 5, 4, 3,
2, 3, 4, 5, 4,
1, 2, 3, 4, 5), ncol = 5)
## square root matrix
Sig.5 <- eigen(Sig)$vectors %*% diag(sqrt(eigen(Sig)$values)) %*% t(eigen(Sig)$vectors)
Sig.5%*%Sig.5
inv.Sig.5 <- solve(Sig.5)
Y <- mvrnorm(n, mu, Sig) # sample mvn data
z <-inv.Sig.5%*%t(sweep(Y, 2, mu))
round(rowMeans(z), 2)
round((cov(t(z))), 2)
## Cholesky Decomposition
T<-chol(Sig)
z <-solve(t(T))%*%t(sweep(Y, 2, mu))
round(rowMeans(z), 2)
round((cov(t(z))), 2)
# increase sample size
Y <- mu, Sig) # sample mvn data
z <-inv.Sig.5%*%t(sweep(Y, 2, mu))
round(rowMeans(z), 2)
round((cov(t(z))), 2)

end video 2

begin video 3

library(MASS) # load library MASS
n <- 1000
mu <- c(10, 5, 7, 9, 20)
Sig <- matrix(c(
5, 4, 3, 2, 1,
4, 5, 4, 3, 2,
3, 4, 5, 4, 3,
2, 3, 4, 5, 4,
1, 2, 3, 4, 5), ncol = 5)

Y <- mvrnorm(n, mu, Sig) # sample mvn data
z.z <- mahalanobis(Y, mu, Sig)

hist(z.z, probability = TRUE, xlab="z'z",
main="Multivariate Normal to Chi-square Distribution")
lines(seq(0, 20, .01), dchisq(seq(0, 20, .01), ncol(Sig)), col='red', lwd=2)
legend("topright", c("Chi-sq Dist w/ df=5"), lty=1, col='red')

library(MASS) # load library MASS
n <- 1000
mu <- c(10, 5, 7, 9, 20)
Sig <- matrix(c(
5, 4, 3, 2, 1,
4, 5, 4, 3, 2,
3, 4, 5, 4, 3,
2, 3, 4, 5, 4,
1, 2, 3, 4, 5), ncol = 5)
Y <- mvrnorm(n, mu, Sig) # sample mvn data

h1 <- hist(Y[, 1], breaks=25, plot=F)
h2 <- hist(Y[, 2], breaks=25, plot=F)
top <- max(h1$counts, h2$counts)
k <- kde2d(Y[, 1], Y[, 2], n=25)
# margins
oldmar <-par()$mar
par(mar=c(3, 3, 1, 1))
layout(matrix(c(2, 0, 1, 3), 2, 2, byrow=T), c(3, 1), c(1, 3))
image(k) #plot the image
par(mar=c(0, 2, 1, 0))
barplot(h1$counts, axes=F, ylim=c(0, top), space=0, col='red')
par(mar=c(2, 0, 0.5, 1))
barplot(h2$counts, axes=F, xlim=c(0, top), space=0, col='red', horiz=T)
par(mar = oldmar)

end video 3

begin video 4

library(MASS) # load library MASS
n <- 1000
mu <- c(1, 2, 0, 3)
Sig <- matrix(c(
4, 0, 1, 3,
0, 4, 1, 1,
1, 1, 3, 1,
3, 1, 1, 9), ncol = 4)

(Sig.yy <- Sig[c(1, 3), c(1, 3)])
(Sig.yx <- Sig[c(1, 3), c(2, 4)])
(Sig.xy <- t(Sig.yx))
(Sig.xx <- Sig[c(2, 4), c(2, 4)])

Cond.Sig <- Sig.yy -
fractions(Cond.Sig)

(mu.y <- mu[c(1, 3)])
(mu.x <- mu[c(2, 4)])

begin video 4