[R] how to use mle with a defined function
Ben Bolker
bolker at ufl.edu
Tue Jul 3 17:10:26 CEST 2007
Lin Pan <linpan1975 <at> yahoo.com> writes:
>
>
> 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'
Can you give a self-contained example (e.g., make up a small
matrix?) There are a few places where I don't understand what
you're doing:
- is the WW above a typo for W?
- if omega is fixed, why calculate its (log) determinant
every time the function is called? (this shouldn't change
the answer, just slow things down)
- is X a single vector or a design matrix?
it seems like you might want lm(U~W-1) instead of lm(U~W)
- it doesn't look like your parameters are used in the
function at all -- are they really parameters for calculating
omega? [because of this, I can get as far as computing
the Hessian, and then I get "system is exactly singular"]
[ off topic: somewhat disturbingly, the CAPTCHA
word for posting from Gmane was mother f***ers ... ]
More information about the R-help
mailing list