# [R] how to use mle with a defined function

Lin Pan linpan1975 at yahoo.com
Tue Jul 10 00:06:13 CEST 2007

```Hi all,

sorry for the misleading in the previous email. here is my function to
calculate the maximum likelihood function for a multinormial distribution:

mymle <- function (sigmaX, sigmaY, constraints, env){

# build omega
omegaX = abs(sigmaX) * kin + abs(env) * diag(1.0, n, n)
omegaY = abs(sigmaY) * kin + abs(env) * diag(1.0, n, n)
omegaXY = (sqrt(sigmaX * sigmaY) - abs(constraints)) * kin
omega = array(0, c(2*n, 2*n))
for(i in 1:n){
for(j in 1:n){
omega[i, j] = omegaX[i, j]
omega[i+n, j+n] = omegaY[i, j]
omega[i+n, j] = omegaXY[i, j]
omega[i, j+n] = omegaXY[i, j]
}
}

# obtain det of omega
odet = unlist(determinant(omega))[1]

# Cholesky decomposition
C = chol(omega)

# beta parameter estimates
newY = t(C)%*%Y
newX = t(C)%*%X

# maximum likelihood estimates
Z = Y - X%*%(lm(newY~newX-1)\$coef)
V = solve(t(C), Z)

# compute the -2log-likelihood
square = t(V)%*%V
0.5*odet + 0.5*square
}

here kin, n, X and Y are known, for example

> kin
[,1]       [,2]       [,3]       [,4]       [,5]       [,6]
[1,] 1.02276611 0.04899597 0.05076599 0.06727600 0.05561066 0.05561066
[2,] 0.04899597 1.02087402 0.11497498 0.07623291 0.06423950 0.06423950
[3,] 0.05076599 0.11497498 1.02291870 0.07189941 0.11162567 0.11162567
[4,] 0.06727600 0.07623291 0.07189941 1.03277588 0.05522156 0.05522156
[5,] 0.05561066 0.06423950 0.11162567 0.05522156 1.03971863 0.54769897
[6,] 0.05561066 0.06423950 0.11162567 0.05522156 0.54769897 1.03971863
> Y
[1] 0.4054651 0.6931472 0.7884574 0.6931472 0.5306283 0.5306283 3.1696856
3.5467397 3.5862929 2.5494452 3.1354942 3.2188758
> X
[,1] [,2] [,3] [,4]
[1,]    1   69    0    0
[2,]    1   65    0    0
[3,]    1   50    0    0
[4,]    1   54    0    0
[5,]    1   48    0    0
[6,]    1   42    0    0
[7,]    0    0    1    1
[8,]    0    0    1    2
[9,]    0    0    1    2
[10,]    0    0    1    2
[11,]    0    0    1    1
[12,]    0    0    1    1
> n
[1] 6

when I call the function
fit <- mle(corr_add, start=list(sigmaX=0.855, sigmaY=0.5, constraints=0.15,
env=0.199),
method = "L-BFGS-B", lower=c(0.001, 0.001, 0.001, 0.001),
upper=c(1000,1000,1000,1000))

it always gave me error message something like
"Error in chol(omega) : the leading minor of order 8 is not positive
definite"

I checked the eigenvalues at each iteration and found that when it stopped
there existed negative eigenvalues. I don't know why this is not working
since omega supposed to be positive definite at each iteration. Any hints
would be very appreciated.

Lin

Lin Pan wrote:
>
> Hi all,
>
> I am trying to use mle() to find a self-defined function. Here is my
> function:
>
> test <- function(a=0.1, b=0.1, c=0.001, e=0.2){
>
> # omega is the known covariance matrix, Y is the response vector, X is the
> explanatory matrix
>
> odet = unlist(determinant(omega))[1]
>
> # do cholesky decomposition
>
> C = chol(omega)
>
> # transform data
>
> U = t(C)%*%Y
> WW=t(C)%*%X
>
> beta = lm(U~W)\$coef
>
> Z=Y-X%*%beta
> V=solve(t(C), Z)
>
> 0.5*odet + 0.5*(t(V)%*%V)
>
> }
>
> and I am trying to call mle() to calculate the maximum likelihood
> estimates for function (0.5*odet+0.5*(t(V)%*%V)) by
>
>
> But I get the following error message:
>
> Error in optim(start, f, method = method, hessian = TRUE, ...) :
>         (list) object cannot be coerced to 'double'
>
> I am pretty sure that the matrices, parameters etc are numerical before
> importing to the function. But why I still get such error message? Could
> anybody give some help on this? thanks a lot.
>
>
> Lin
>

--
View this message in context: http://www.nabble.com/how-to-use-mle-with-a-defined-function-tf4015002.html#a11511446
Sent from the R help mailing list archive at Nabble.com.

```