STA313 F 2004 Handout 3

Numerical MLE for the Multivariate normal with mvnLL

First, illustrate some commands.

> Ymat[1:10,] # First 10 rows of Ymat
               y1         y2
 [1,]  0.46950842  0.6999623
 [2,]  0.75046599  0.3073191
 [3,] -0.03310517  1.6231372
 [4,]  1.95780725  3.3519020
 [5,]  0.68899300  5.1409713
 [6,]  2.95751217  2.9669962
 [7,]  1.79292404 -0.5733391
 [8,] -0.51203666 -0.4667499
 [9,] -1.40452153 -1.9487259
[10,] -2.41753657  0.1730442

> dim(Ymat)
[1] 200   2

> Sigma
          y1        y2
y1 1.8429904 0.8355763
y2 0.8355763 1.8324584

>
> siginv <- solve(Sigma) ; siginv
          [,1]       [,2]
[1,]  0.6840044 -0.3118968
[2,] -0.3118968  0.6879357

> Ymat[1,]
       y1        y2
0.4695084 0.6999623
>
> Ymat[1,]%*%siginv
          [,1]      [,2]
[1,] 0.1028298 0.3350909
>
> Ymat[1,]%*%siginv%*%cbind(Ymat[1,])

            [,1]
  [1,] 0.2828305
>
> mean(Ymat[,1]) # Mean of col 1
[1] 1.096919
>
> apply(Ymat,2,mean)
      y1       y2
1.096919 2.069276
>

> nn <- dim(Ymat)[1] ; nn
[1] 200

Here is the function definition. It's in the file mvnLL.R.

# mvnLL.R : Defines the function mvLL
#           source("mvnLL.R"), probably in a program defining another function.
#
mvnLL <- function(Sigma,Ymat,check=T,centered=F)
#        -2 LL for a multivariate normal: Function of the covariance matrix
    {
    p <- dim(Ymat)[2]
    if(check)
        {
        p1 <- dim(Sigma)[1] ; p2 <- dim(Sigma)[2]
        if(p1 != p2)
            {cat("     Sigma = \n")
            print(Sigma)
            stop("Sigma must be square!")
            }
        if(p1 != p) stop("Num cols of Y must = num cols of Sigma.")
        } # End checking of input
    deter <- det(Sigma)
    # If determinant is neg or very near zero, return +Infinity
    if(deter < 4e-16) mvnLL <- Inf else
       {
       siginv <- solve(Sigma)
       nn <- dim(Ymat)[1]
       const <- nn*p*log(2*3.1415926)
       sumq <- 0 # Initializing sum of quadratic forms.
       for(i in 1:nn) sumq <- sumq +
            as.numeric( Ymat[i,]%*%siginv%*%cbind(Ymat[i,]) )
       if(! centered)
           {
           ybar <- cbind(apply(Ymat,2,mean)) # Y-bar is a column vector
           sumq <- sumq - nn*t(ybar)%*%siginv%*%ybar
           }
       mvnLL <- as.numeric( nn*log(deter) + const + sumq ) # Make it a number
       } # End if the determinant is okay
    mvnLL
    } # End of function mvnLL

Need to wrap mvnLL up in another function to minimize.

# simplecov.R
simplecov <- function(theta,datta) # For a 2x2
    {
    sig11 <- theta[1] ; sig12 <- theta[2] ; sig22 <- theta[3]
    covmat <- rbind( c(sig11,sig12),
                     c(sig12,sig22)   )
    simplecov <- mvnLL(covmat,datta)
    simplecov
    } # End of simplecov

Test the functions.

> #                    Illustrate MLE with mvnLL
> # Simulate some data
> n <- 200 ; set.seed(9999)
> f1 <- rnorm(n,1,1) ; f2 <- rnorm(n,2,1) ; e <- rnorm(n)
> x1 <- f1+e ; x2 <- f2+e ; XX <- cbind(x1,x2)
> # Define the functions
> source("mvnLL.R") ; source("simplecov.R")
> out1 <- nlm(simplecov,c(1,0,1),datta=XX)[1:2] ; out1
Warning messages:
1: NA/Inf replaced by maximum positive value
2: NA/Inf replaced by maximum positive value
3: NA/Inf replaced by maximum positive value
$minimum
[1] 1330.040

$estimate
[1] 2.003452 1.088755 1.914238

> sig11 <- out1$estimate[1] ; sig12 <- out1$estimate[2]
> sig22 <- out1$estimate[3]
> Sigmahat <- rbind(c(sig11,sig12),
+                   c(sig12,sig22))
> Sigmahat # Numerical MLE
         [,1]     [,2]
[1,] 2.003452 1.088755
[2,] 1.088755 1.914238
> Sighat <- (n-1)/n * var(XX)
> Sighat   # Theoretical MLE
         x1       x2
x1 2.003457 1.088759
x2 1.088759 1.914241
> #
> out1$minimum # Minimum value of -2LL from numerical search
[1] 1330.040
> # Now evaluate (simplified) -2LL at the the Theoretical MLE
> p <- 2 ; n*p*(1+log(2*3.14159)) + n*log(det(Sighat))
[1] 1330.04
> #
> # Now try some different starting values
> #
> testme1 <- nlm(simplecov,c(2,1,2),datta=XX)[1:2] ; testme1
Warning messages:
1: NA/Inf replaced by maximum positive value
2: NA/Inf replaced by maximum positive value
$minimum
[1] 1330.040

$estimate
[1] 2.003469 1.088770 1.914252


> testme2 <- nlm(simplecov,c(6,0,23),datta=XX)[1:2] ; testme2
There were 14 warnings (use warnings() to see them)
$minimum
[1] 1330.040

$estimate
[1] 2.003453 1.088756 1.914237

> testme3 <- nlm(simplecov,c(-1,0,3),datta=XX)[1:2] ; testme3
Warning messages:
1: NA/Inf replaced by maximum positive value
2: NA/Inf replaced by maximum positive value
3: NA/Inf replaced by maximum positive value
4: NA/Inf replaced by maximum positive value
$minimum
[1] 1.797693e+308

$estimate
[1] -1  0  3

> #
> # Starting at the identity matrix seems like a good practice for now.
> #
>