[R] how to use mle with a defined function

[Lin Pan]

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
> 
> result = mle(test, method="Nelder-Mead")
> 
> 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.