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.
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
