pois_reg_sarah.qmd

---
title: "pois_reg"
format: html
---

The clean versions of the score and information functions unfortunately do not live in matrix syntax.


```{r}
#score function
U <- function(beta, x, y) {
  
  #initialize lambda, score, sample size
  l <- exp(beta[1] + beta[2]*x)
  score <- rep(0, times = 2)
  
  #fill first element of score
  score[1] <- sum(y - l)
  
  #fill second element of score
  score[2] <- sum(x * (y - l))
  
  #return score
  score
}

#Fisher information
info <- function(beta, x, y) {
  
  #initialize Hessian, lambda, sample size
  Hessian <- matrix(data = c(0,0,0,0), nrow = 2, ncol = 2)
  l <- exp(beta[1] + beta[2]*x)
  
  #fill top left entry in Hessian
  Hessian[1,1] <- - sum(l)
  
  #fill mixed entries in Hessian
  Hessian[1,2] <- Hessian[2,1] <- - sum(x * l)
  
  #fill in bottom right entry in Hessian
  Hessian[2,2] <- - sum(x ^ 2 * l)
  
  #return Hessian
  Hessian
}
```

```{r}
find_beta <- function(start_guess, x, y, maxtol = 1E-4, maxiter = 1E3, verbose = FALSE) {
  beta_curr <- start_guess
  diff <- 1
  iterations <- 0
  while(diff > maxtol & iterations <= maxiter){
    beta_next <- beta_curr - solve(info(beta_curr, x, y)) %*% U(beta_curr, x, y)
    diff <- max(abs(beta_next[1] - beta_curr[1]), abs(beta_next[2] - beta_curr[2]))
    iterations <- iterations + 1
    beta_curr <- beta_next
    if (verbose) {
      print(paste("beta_next:", beta_next))
      print(paste("diff:", diff))
      print(paste("iterations:", iterations))
    }
  }
  
  if (diff > maxtol & iterations >= maxiter) {
    conv_msg = "Hit maxiter without convergence."
  } else {
    conv_msg = "Yay! We converged."
  }
  
  return(list(estimates = beta_curr, 
              convergence = conv_msg))
}
```

Now, to test.
```{r}
set.seed(1031)

x <- c(1,2,3,4)
beta <- c(5,6)
lambda <- exp(5 + 6 * x)
y <- rpois(4, lambda)

find_beta(c(4.9, 5.9), x, y) #problem

summary(glm(y ~ x, family = "poisson")) #did not converge, maybe that's why?

```

```{r}
b0 = .5
b1 = .25
x = runif(100, min = 0, max = 1)
lambda = exp(b0 + b1*x)
y = rpois(100, lambda = lambda) 

set.seed(1031) ## happy birthday!
find_beta(c(0, 0), x, y)

summary(glm(y ~ x, family = "poisson")) #did not converge, maybe that's why?
```