[R] Log-Likelihood 3d-plot and contourplot / optim() starting values
voodooochild@gmx.de
voodooochild at gmx.de
Wed Jan 25 18:04:05 CET 2006
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
More information about the R-help
mailing list