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))