Find the bug in the following code that seeks to use NR method to optimize alpha and beta for a beta distribution: #Prep
set.seed(0820)
library(tidyverse)
cleaned_dataset <- read_csv("cleaned_dataset.csv")

cleaned_dataset_sub <- cleaned_dataset %>% 
  filter(G >= 41, `FT%` != 0, `FT%` != 1) #Note throw out 5 players at extremes for errors

#Scrape test FT%
obs <- cleaned_dataset_sub %>% 
  filter(season <= 2013) %>% 
  select(`FT%`) %>% 
  unlist() %>% 
  unname()

test <- cleaned_dataset_sub %>% 
  filter(season > 2013) %>% 
  select(`FT%`) %>% 
  unlist() %>% 
  unname()

#log gamma function
log_gamma <- function(n) { #converted to log scale since numbers so big
  log(factorial(n-1))
}

#log beta function
log_beta_func <- function(a, b) {
  log_gamma(a) + log_gamma(b) - log_gamma(a+b)
}

#beta density values
beta_d <- function(x, a, b) {
  exp((a-1) * log(x) + (b-1) * log(1-x) - log_beta_func(a,b))
}

#minus log likelihood for beta distribution
minus_log_likelihood <- function(par, data) {
  alpha <- par[1] 
  beta <- par[2]
  -sum(log(beta_d(data, alpha, beta)))
}

#First derivative of the log likelihood
ll_gradient <- function(par, data) {
	alpha <- par[1]
	beta <- par[2]
	dalpha <- sum(1/beta(alpha, beta) * (log(data) + (alpha-1)/data - (beta-1)/(1-data)))
  dbeta <- sum(1/ beta(alpha, beta) * ((alpha-1)/data - (beta-1)/(1-data)))
  
	matrix(c(dalpha,
	       dbeta), ncol = 1)
}



#Second derivative of the log likelihood
ll_hessian <- function(par, data) {
	alpha <- par[1]
	beta <- par[2]
  b <- beta(alpha, beta)
  dalpha2 <- -sum(1/b * (log(data) + (beta-1)/(1-data) - alpha/data)^2 - 1/b/data^2)
  dbeta2 <- -sum(1/b * (log(1-data) + (alpha-1)/data - beta/(1-data))^2 - 1/b/(1-data)^2)
  dalpha_dbeta <- sum(1/b * ((1/(1-data) + 1/data) * (log(1-data) + (alpha-1)/data - (beta-1)/(1-data))))
  H <- matrix(c(dalpha2, dalpha_dbeta, dalpha_dbeta, dbeta2), ncol = 2)
  H
}

#Choose arbitrary starting values to inital function
#Newton-Raphson
a_0 <- 8
b_0 <- 2
theta <- c(a_0, b_0)

l <- -minus_log_likelihood(theta, obs)
g_prime <- ll_gradient(theta, obs)
g_prime_prime <- ll_hessian(theta, obs)
eps = 0.00000001 #initial acceptable epsilon value

while(sum(abs(g_prime))>eps) #convergence
{
  theta = theta - solve(g_prime_prime, g_prime)
  #multivariate optimization technique
  l <- -minus_log_likelihood(theta, obs)
  g_prime <- ll_gradient(theta, obs)
  g_prime_prime <- ll_hessian(theta, obs)
}
# 
c(alpha_mle = theta[1], beta0_mle = theta[2])