##  This is a practicaly demonstration of how to obtain unbiased estimates from
##  a Poisson regression in the presence of random noise


##  The outcome variable y is drawn from a Poisson distribution with lambda
##  equal to exp(1 + 2 * x). Where x is a random draw from the standard normal
##  distribution.

##  Data simulation
n <- 10000 # number of draws
x <- rnorm(n)
y <- rpois(n, lambda = exp(1 + 2 * x))

##  Estimating the parameters using the default Poisson regression model in R
glm(y ~ x , family = poisson) ## we obtain correct estimates 

##  Adding random noise: z is the same as y but with an amount of random noise
##  u added. u is drawn from the Poisson distribution with lambda = 3
u <- rpois(n, lambda = 3)
z <- y + u

##  Estimating again using the default in R
glm(z ~ x , family = poisson) ## our parameters are no longer unbiased

##  The correct routine:
##  First we specify the log-likelihood function

pois.loglik <- function(beta){
  lambda <- exp(cbind(1, x) %*% beta) + 3 # This is the correct expected lambda
  log.lik <- sum(dpois(z, lambda, log = T)) # Obtain the log likelihood
  return(log.lik)
}

##  Then we find the values of beta that maximise the log-likelihood using a 
##  optimiser function
optim(par = 3:2, #initial guesses for beta
      fn = pois.loglik, #function that produces the log-likelihood 
      control = list(fnscale = -1) #final option to say we wish to maximise 
      )
##  The outcome under $par are now the correct values of beta.