R optimization with optim

Question

I have no idea what went wrong in the following code, after suffering for a long time,possibly whole day... i really hope someone can find the problem...Thanks a lot in advance ....I don't know why R is giving me the same initial value

lambda<-1 

likelihood<- function(lambda) {
  T6<- exp(lambda*2.82)/(exp(lambda*2.82)+exp(lambda*0.28))
  T5<- exp(lambda*2.69)/(exp(lambda*2.69)+exp(lambda*((T6*0.38)+((1-T6)*2.92))))
  T4<- exp(lambda*2.52)/(exp(lambda*2.52)+exp(lambda*((T5*0.51)+((1-T5)*2.82))))
  T3<- exp(lambda*2.30)/(exp(lambda*2.30)+exp(lambda*((T4*0.68)+((1-T4)*2.69))))
  T2<- exp(lambda*2.00)/(exp(lambda*2.00)+exp(lambda*((T3*0.90)+((1-T3)*2.52))))
  T1<- exp(lambda*1.60)/(exp(lambda*1.60)+exp(lambda*((T2*1.20)+((1-T2)*2.30))))

  l6<-(T6^0)*((1-T6)^281)
  l5<-(T5^0)*((1-T5)^281) 
  l4<-(T4^2)*((1-T4)^(281-2)) 
  l3<-(T3^19)*((1-T3)^(281-19))
  l2<-(T2^93) *((1-T2)^(281-93))
  l1<-(T1^166)*((1-T1)^(281-166)) 

         -(l1*l2*l3*l4*l5*l6)
         }

optim(lambda,likelihood)$par

LyzandeR LyzandeR · Accepted Answer · 2015-05-17T16:16:14

There is nothing wrong with the code (in terms of a programming error) . If you feel that it shouldn't return the upper limit then you probably made a mistake with the definition of the function in your code.

First of all when you only have one parameter to estimate your method should be Brent like this:

> optim(lambda,likelihood, method='Brent', lower=1.5, upper=200)$par
[1] 200

The reason it always shows the upper limit as the correct one is because your function is only defined between -100 - 300 (approximately) and its value is always zero as you can see below (irrespective of the value of the lambda). When this happens optim always shows the upper limit as the right value. You can see your function plotted below:

plot(likelihood, from=-100, to=1000, xlab='lambda')

enter image description here

R optimization with optim

1 Answers