#####  https://cmdlinetips.com/2019/03/introduction-to-maximum-likelihood-estimation-in-r/
# modified by Pedro Quintana-Ascencio and Federico Lopez Borghesi for Methods II class 2023

data <- rpois(n=1000, lambda =5 )
df_Poiss <- data.frame(data)

###############################################################
# plotting  an approximation of the population and a sample
###############################################################
par(mfrow=c(2,2))
b <- seq(0,20,0.5)
data <- rpois(n=20000, lambda =5 )
df_Poiss <- data.frame(data)
hist(df_Poiss$data,breaks = b,main="stat pop (aprox")

################################
## SAMPLED DATA

datap <- rpois(n=10, lambda =5 )
################################

df_Poiss <- data.frame(datap)
hist(df_Poiss$datap,breaks = b,main="sample")


###############################################################
#likelihood for a single data point with a single parameter
###############################################################

data_chosen <- datap[3]
likelihood <- dpois(data_chosen,lambda=5)

###############################################################
#likelihood for a single data point with multiple parameters 
###############################################################

#par(mfrow=c(1,2))
likelihood <- dpois(data_chosen,lambda=seq(1,20,1))
plot(likelihood, pch=16,bty="n",main= paste("for a single data",data_chosen))
abline(v=data_chosen)
log_likelihood <- dpois(data_chosen,lambda=seq(1,20,1),log=TRUE)
plot(log_likelihood, pch=16,bty="n",main= paste("for a single data",data_chosen))
abline(v=data_chosen)
abline(v=data_chosen)

###############################################################
#likelihood for multiple data points with a multiple parameter
# (lambdas)
###############################################################

ll_Poiss <- function(lambda,y){
  llh <- sum(dpois(y,lambda,log=TRUE))
  return(llh)
}

lambdas <- seq(1,15,by=0.5) # parameters

ll <- sapply(lambdas,function(x){ll_Poiss(x,data)})

par(mfrow=c(1,1))
df <- as.data.frame( t(rbind(ll,lambdas)) )
plot(df$lambdas,df$ll,pch=16,type="l", main ="log_likelihood for multiple lambdas for sampled data",
     bty="n",xlim=c(0,18))
abline(v=5) #lambda used to generate the data


#############################################################
# Estimating parameter "PROBABILITY" using STAN
#############################################################

library(rethinking)
library(rstan)
model1 <- map2stan(
  alist(
    datap ~ dpois(lambda),
    log(lambda) <- t,
    t ~ dnorm(0,20)
  ),
  data = df_Poiss,chains =3,
  iter=6000,warmup=2000
)

precis(model1,digits=1,depth=1)
post <- extract.samples(model1)
plot(density(exp(post$t)),type="l", main="probability of lambda for the sampled data",
     bty="n",xlim=c(0,18))
abline(v=5) #lambda used to generate the data