[R] Constraint maximum (likelihood) using nlm

[Benjamin Dickgiesser]

Nevermind, I had a small mistake in the negll function.

Thank you anyways.

On 2/17/07, Benjamin Dickgiesser <dickgiesser at gmail.com> wrote:
> Hi,
>
> I'm trying to find the maximum (likelihood) of a function. Therefore,
> I'm trying to minimize the negative likelihood function:
>
> # params: vector containing values of mu and sigma
> # params[1] - mu, params[2]- sigma
> # dat: matrix of data pairs y_i and s_i
> # dat[,1] - column of y_i , dat[,2] column of s_i
> negll <- function(params,dat,constant=0)
> {
>         for(i in 1:length(dat[,1]))
>         {
>                 llsum <- log( params[2]^2 + dat[i,2]^2) +
>                 (( dat[i,1] - params[1])^2/ (params[2]^2 + dat[i,2]^2))
>         }
>         ll <- -0.5 * llsum  + constant
>         return(-ll)
> }
>
> Using (find data attached):
> data.osl <- read.table("osl.dat",header=TRUE)
> data.matrix <- as.matrix(data.frame(data.osl$de,data.osl$se))
> nlm(negll,c(0.75,0.5),dat=data.matrix,iterlim=200)
>
> I get estimates for mu and sigma of:
>  3.629998e+00 -4.975368e-07
>
> However, sigma obviously has to be >= 0.
>
> Therefore I am trying to transform sigma:
>
> negll.trans <- function(params,...)
> {
>         params[2] <- log(params[2])
>         negll(params,...)
>
> }
>
> where
> nlm(negll.trans,c(0.75,0.5),dat=data.matrix)
> give me the estimates
> 3.63 1.00
> i.e. mu = 3.63 and sigma = exp(1)
>
> I am not confident that this is the correct answer looking at the graphs:
> par(mfrow=c(2,1))
> plot(osl$de,osl$se,xlab="equivalent dose",ylab="standard errors",
>      main="Equivalent Dose vs. Standard Errors",xlim=c(0,4))
> hist(osl$de,xlab="equivalent dose",
>      main="Histogram of Equivalent Dose",xlim=c(0,4))
>
> I want to stick with using nlm to minimize the function but can't find
> an error in what I am doing.
>
> I'd appreciate your help!
>
> Thank you
> Ben
>
>