voodooochild@gmx.de
2006-Jan-25 17:04 UTC
[R] Log-Likelihood 3d-plot and contourplot / optim() starting values
Hello, i have coded the following loglikelihood-function # Log-Likelihood-Funktion loglik_jm<-function(N,phi,t) { n<-length(t) i<-seq(along=t) s1<-sum(log(N-(i-1))) s2<-phi*sum((N-(i-1))*t[i]) n*log(phi)+s1-s2 } # the data t<-c(7,11,8,10,15,22,20,25,28,35) # now i want to do a 3d-plot and a contourplot in order to see at which values of N and phi the loglikelihood function becomes zero. # i do this in order to get an idea where the starting values for a optimization for the mle of N could be # 3dplot and contourplot phi<-seq(0,1,length=50) N<-seq(101,110,length=50) z<-outer(N,phi,loglik_jm(N,phi,t)) persp(phi,N,z, theta = 30, phi = 30, expand = 0.5, col = "lightblue") contourplot(z~N*phi) # but i get some error messages, i don't know why? # if you are interested, the mle function for N is ll2 <- function(N,t) { i<-seq(along=t) n<-length(t) s1<-sum(1/(N-(i-1))) s2<-n/(N-(1/sum(t[i]))*(sum((i-1)*t[i]))) (s1-s2) } # you get this function as usual, if you set the loglikelihood equal zero and differentiate for N # i take the squares of ll2 in order to get the minimum easier ll3<-function(N,t) { ll2(N,t)^2 } # then i do an iteration to get the mle of N xmin<-optim(10,ll3,method="BFGS",control=list(reltol=.Machine$double.eps),t=t) # the problem is that this method only works well, if the starting values for optim() are very close to the real value! # so i got the idea with contourplot() and persp() to see where good starting values could be. i would be very thankful if anyone could give me some advice! regards Andreas