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
>
>