## generate random data
y <- rpois(25, 2)
x <- rnorm(25) + (log(replace(y, y == 0, 0.25)))

mod1 <- optim(c(1, 1), function(p)
    with(list(mu=p[1] + p[2] * x), -sum(y * mu - exp(mu))))

mod2 <- glm(y ~ x, family='poisson', data=data.frame(y=y, x=x))

## check that the parameters coincide
mod1$par
coefficients(mod2)

## eta = a + bx
eta <- cbind(1, x) %*% mod1$par
tail(eta)
tail(predict(mod2))

## yhat = exp(eta)
yhat <- exp(eta)
tail(yhat)
tail(fitted.values(mod2))

## log-likelihood
sum(dpois(y, yhat, log=TRUE))
logLik(mod2)[1]

## deviance: D = 2(logLik_fullmodel - logLik_model)
2 * (sum(dpois(y, y, log=TRUE)) - sum(dpois(y, yhat, log=TRUE)))
mod2$deviance

tail(2 * (yhat - y - y * (eta - log(y))))
tail(resid(mod2)^2)
sum(resid(mod2)^2)

barplot(rbind(dpois(0:4, y[1], log=TRUE),
              dpois(0:4, yhat[1], log=TRUE)),
        beside=TRUE, names.arg=paste0('y=', 0:4),
        legend.text=c('saturated', 'model'), args.legend=list(x='bottom'))