set.seed(322)

n<-100
X<-runif(n,-10,10)
b0<-1.0
b1<-0.5
mu<-exp(b0+b1*X)
Y<-rpois(n,mu)

glm.p<-glm(Y~X,family=poisson)

#Initialization
Xmat<-cbind(rep(1,n),X)

beta.val<-c(0,0)

beta.val<-c(log(mean(Y)),0)


eta<-Xmat %*% beta.val
muvec<-exp(eta)
W<-diag(as.vector(muvec))
Z<-eta+(Y-muvec)/muvec

eqns<-t(Xmat) %*% (Y-muvec)

#Recursion
istep<-0

while(sum(eqns^2) > 1e-8){

	beta.val<-solve(t(Xmat) %*% W %*% Xmat) %*% t(Xmat) %*% W %*% Z
	eta<-Xmat %*% beta.val
	muvec<-exp(eta)
	W<-diag(as.vector(muvec))
	Z<-eta+(Y-muvec)/muvec
	eqns<-t(Xmat) %*% (Y-muvec)	
	istep<-istep+1
	print(c(istep,sum(eqns^2)))
}

print(istep)	
print(t(beta.val))
print(sqrt(diag(solve(t(Xmat) %*% W %*% Xmat))))

#Compare results
summary(glm.p)