From: <voodooochild_at_gmx.de>

Date: Thu 26 Jan 2006 - 04:04:05 EST

persp(phi,N,z, theta = 30, phi = 30, expand = 0.5, col = "lightblue") contourplot(z~N*phi)

}

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Jan 26 04:41:08 2006

Date: Thu 26 Jan 2006 - 04:04:05 EST

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

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Jan 26 04:41:08 2006

*
This archive was generated by hypermail 2.1.8
: Thu 26 Jan 2006 - 08:09:38 EST
*