rm(list = ls(all = TRUE))
graphics.off()

## --------------------------------- 多元正态分布的基本性质 ---------------------------------
x = 0
mu = -0.8 * x
sigma = sqrt(1.36)
curve(dnorm(x, mu, sigma), -10, 5, col=1, xlab='', ylab='', 
      lwd=2, ylim=c(0, 0.35), 
      main='Conditional Normal Densities f(X2|X1)')
for (i in 1:6) {
  x = x + 0.85
  mu = -0.8 * x
  curve(dnorm(x, mu, sigma), -10, 5, col=i+1, xlab='', ylab='', 
        add=TRUE, lty=i+1, lwd=2)
}

## Conditional Approximations
x <- c(230, 125, 200, 109, 181, 99, 55, 107, 165, 97, 105, 98, 
       150, 115, 85, 71, 97, 120, 0, 82,
       192, 100, 150, 103, 181, 80, 85, 111, 189, 90, 120, 93, 
       172, 95, 110, 86, 170, 125, 130, 78)
pullover <- matrix(x, ncol = 4, byrow = TRUE)
pullover <- as.data.frame(pullover)
colnames(pullover) <- c("X1", "X2", "X3", "X4")
pullover
mu = apply(pullover, 2, mean)
mu
Sigma = var(pullover) * (9 / 10)
Sigma

mu_2 = as.matrix(mu[2:4])
Sigma_21 = as.matrix(Sigma[1, 2:4])
Sigma_22 = as.matrix(Sigma[2:4, 2:4])
t(Sigma_21) %*% solve(Sigma_22)
mu[1] - t(Sigma_21) %*% solve(Sigma_22) %*% mu_2

Sigma[1, 1] - t(Sigma_21) %*% solve(Sigma_22) %*% Sigma_21

t(Sigma_21) %*% solve(Sigma_22) %*% Sigma_21 / Sigma[1, 1]

cor(pullover)

mu_1 = as.matrix(mu[1:2])
mu_2 = as.matrix(mu[3:4])
Sigma_11 = as.matrix(Sigma[1:2, 1:2])
Sigma_12 = as.matrix(Sigma[1:2, 3:4])
Sigma_21 = as.matrix(Sigma[3:4, 1:2])
Sigma_22 = as.matrix(Sigma[3:4, 3:4])
Sigma_12 %*% solve(Sigma_22)
mu_1 - Sigma_12 %*% solve(Sigma_22) %*% mu_2

Sigma_11 - Sigma_12 %*% solve(Sigma_22) %*% Sigma_21


## --------------------------------- 球形分布与椭圆分布 ---------------------------------
n = 80
mu1 = 0
mu2 = 0
s1 = 1
s2 = 1
rho = 0.6
x = seq(-3, 3, length = n) * s1
y = seq(-3, 3, length = n) * s2
f = function(x,y){
  (2 * pi * s1 * s2 * sqrt(1-rho^2))^-1 * exp(-0.5 * (1 - rho^2)^-1 * 
                ((x-mu1)^2/s1^2 - 2 * rho * (x - mu1) * ( y - mu2) / (s1 * s2) 
                 + (y - mu2)^2/s2^2))
}
z = outer(x, y, f)

library(ggplot2)
library(plotly)
fig_1 = plot_ly() %>%
  add_surface(x = ~x, y = ~y, z = ~z)
fig_1

## persp(x, y, z, theta = -45, phi = 45)

fig_2 = plot_ly(x = ~x, y = ~y, z = ~z) %>%
  add_contour()
fig_2