STAT 3690 Lecture 08
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit
STAT 3690 Lecture 08
2022
Assumptions
● Model: x1 , . . . , xn MVNp (μ, 3), n > p
● Parameter space: Θ = {(μ, 3) l μ e Rp , 3 e Rp×p , 3 > 0}
Method of moments (MM) estimators for (μ , 、)
● No requirement on normality
● Steps
1. Equate raw moments to their sample counterparts:
i(T) 兮 i(T)
2. Solve the above equations w.r.t. μ and 3 and obtain estimators
i(T) _ T = n_ 1 (n _ 1)S,
where S = (n _ 1)_ 1 (xi _ )(xi _ )T
Maximum likelihood (ML) estimation for parameters of MVN (J&W Sec 4.3)
● Likelihood function
┌ exp - _ (xi _ μ)T 3_ 1 (xi _ μ) 入┐
= ,(2π)np {(1)det(3)}n exp - (xi _ μ)T 3_ 1 (xi _ μ) 入
● Log likelihood
e(μ, 3) = ln{L(μ, 3)} = _ ln(2π) _ ln{det(3)} _ n (xi _ μ)T 3_ 1 (xi _ μ)
i=1
options(digits = 4)
set.seed(1)
n = 1e3L
Mu = matrix(c (3 , 6), ncol = 1 , nrow = 2)
Sigma = matrix(c (10 , 2 ,2 , 5), ncol = 2 , nrow = 2)
X = MASS::mvrnorm(n, Mu, Sigma)
loglik <- 达{≥′t↓μ≥(Mu, Sigma, data = X) {
n = nrow (data)
p = length (Mu)
X_center = sweep(X, 2, Mu)
return(
(-n*p/2)*log(2*pi)+
(-n/2)*log(det(Sigma))+
(- 1/2)*sum (diag(X_center %*% solve(Sigma) %*% t(X_center)))
)
}
grid_xy <- expand.grid(
seq(Mu[1]-2*Sigma[1 ,1] ˆ .5, Mu[1]+2*Sigma[1 ,1] ˆ .5 , length.out = 32),
seq(Mu[2]-2*Sigma[2 ,2] ˆ .5, Mu[2]+2*Sigma[2 ,2] ˆ .5 , length.out = 32))
contours <- purrr::map_df(
seq_len(nrow (grid_xy)),
达{≥′t↓μ≥ (i) {
# Where we will evaluate loglik
mu_obs <- as.numeric(grid_xy[i,])
# Evaluate at the pop covariance
z <- loglik(mu_obs, Sigma, X)
# Output data.frame
data.frame(x = mu_obs[ 1],
y = mu_obs[2],
z = z)
})
# Contour plot
library(tidyverse)
library(ggrepel)
data_means <- data.frame(
x = c (Mu[1], mean (X[,1])),
y = c (Mu[2], mean (X[,2])),
label = c ( "Mu" , "Sample Mean"))
contours %>%
ggplot(aes (x, y)) +
geom_contour(aes (z = z)) +
geom_point(data = data_means, aes (color = label)) +
geom_label_repel(data = data_means, aes (label = label))
# 3d scatter plot
library(scatterplot3d)
with(contours, scatterplot3d(x, y, z))
2022-03-17