#log-likelihood function for poisson distribution

like<-function(lam,n,x){
    -(n)*lam+log(lam)*sum(x)-sum(log(x))-n*log(1-exp(-lam))
}


dlike<- function(lam,n,x){
    -(n)+(sum(x)/lam)-n*exp(-lam)/(1-exp(-lam))
}

dlike2<-function(lam,n,x) {
    -sum(x)/lam^2+n*exp(-lam)/(1-exp(-lam))^2
}

score<-function(lam,n) {
    - n/(lam*(1-exp(-lam)))+(n*exp(-lam)/(1-exp(-lam))^2)
}


modified.newton<- function(fun, derf, derf2, x0,n,x,eps){
    k<- 0
    repeat {
	k<- k + 1
	x1<- x0 - derf(x0,n,x)/derf2(x0,n,x)
	if(abs(x0 - x1) < eps)
	    break
	x0<- x1
	cat("Iteration No: ", k, "Current iterate = ",x1 ,"Current Likelihood = ",fun(x0,n,x),"dl/dx = ",derf(x0,n,x),"dl2/dx2=",derf2(x0,n,x),fill=T)
    }
    x1
}
   
 newton.score<- function(fun, derf, score, x0,n,x,eps){
    k<- 0
    repeat {
	k<- k + 1
	x1<- x0 - derf(x0,n,x)/score(x0,n)
	if(abs(x0 - x1) < eps)
	    break
	x0<- x1
	cat("Iteration No: ", k, "Current iterate = ",x1 ,"Current Likelihood = ",fun(x0,n,x),"dl/dx = ",derf(x0,n,x),"Score=",score(x0,n),fill=T)
    }
    x1
}


# produce a data vector x
x<-rep(1:6,c(1486,694,195,37,10,1))
n<-length(x)
table(x)

# plot the log-likelihood function
ltemp<- seq(0.1,10,0.1)
lhood<- like(ltemp,n,x)
plot(ltemp,lhood,type="l",xlab="parameter values",
	ylab="log-likelihood", main="Log-likelihood function of the truncated Poisson distribution")



# initial value
est<- mean(x)
final.est<-modified.newton(like,dlike,dlike2,est,n,x,0.000001)
final.est
est<- mean(x)
final.score<-newton.score(like, dlike, score, est,n,x,0.000001)
final.score