#Ejemplo del curso 2023
#Trabajan con la densidad
p <- function(x){
  y <- 0.5*(1+theta*x)
  return(y)
}

#Genero un vector de 7 datos
vec=c(-0.47,0.89,0.26,-0.59,0.37,0.54,0.87)

L <- function(theta){ #likelihood
  prod(0.5*(1+theta*vec))
}
Lv <- Vectorize(L)#para dibujar
nlogL <- function(theta){ #loglikelihood negativo porque optim minimiza
  -sum(log(0.5*(1+theta*vec)))
}

nlogLv <- Vectorize(nlogL)#para dibujar
layout(matrix(1:2,ncol = 2))
M=optimize(f=nlogL,interval=c(-1,1))$minimum #estimador MLE
layout(matrix(1:2,ncol = 2))
curve(Lv,from=-1,to=1,n = 1024,xlab = expression(theta),ylab = expression(L(theta)),main="Verosimilitud",lwd=2)
abline(v=M,lty=2,lwd=2,col=2)
legend(x="topleft",cex=0.7,horiz=F, c('Verosimilitud','MLE'),col=c(1:2),lty=c(1,2),lwd=c(1,2))
curve(nlogLv,from = -1,to = 1,n=1024,xlab = expression(theta),ylab = expression(-lnL(theta)),main="-Log Verosimilitud",lwd=2)
abline(v=M,lty=2,lwd=2,col=2)
legend(x="topleft",cex=0.7,horiz=F, c('-Log Verosimilitud','MLE'),col=c(1:2),lty=c(1,2),lwd=c(1,2))


########################
#Descenso por gradiente#
########################
theta0 <- -0.8
alpha <- 0.1
vec=c(-0.47,0.89,0.26,-0.59,0.37,0.54,0.87)

layout(1)
curve(nlogLv,from = -1,to = 1,n=1024,xlab = expression(theta),ylab = expression(-lnL(theta)),main="-Log Verosimilitud",lwd=2)
points(theta0,nlogL(theta0),pch=19,col="red")
text(theta0,nlogL(theta0)+2,"0",col="red")

#Calculamos el gradiente
gradiente <- function(theta){
  -sum(vec/(1+theta*vec))
}

#Iteramos
theta1 <- theta0-alpha*gradiente(theta0)
points(theta1,nlogL(theta1),pch=19,col="blue")
text(theta1,nlogL(theta1)+2,"1",col="blue",cex=0.5)

theta2 <- theta1-alpha*gradiente(theta1)
points(theta2,nlogL(theta2),pch=19,col="blue")
text(theta2,nlogL(theta2)+2,"2",col="blue",cex=0.5)

thetak <- theta2
for(i in 3:100){
  thetakk <- thetak-alpha*gradiente(thetak)
  points(thetakk,nlogL(thetakk),pch=19,col="blue")
  text(thetakk,nlogL(thetakk)+2,i,col="blue",cex=0.5)
  i <- i+1
  thetak <- thetakk
  readline(prompt="Press [enter] to continue")
}

#Genero varios vectores con n=7 datos
#theta=0.5 el verdadero parametro
theta <- 0.5
#Definimos la densidad p(x)
p <- function(x){
  y <- 0.5*(1+theta*x)
  return(y)
}
#Para usar algunas funciones es necesario vectorizarla
pv <- Vectorize(p)
curve(pv,from = -1,to = 1,xlab = "x",ylab = "p(x)",lwd=2,main="Densidad")
M <- optimize(p, interval=c(-1, 1), maximum=TRUE)$objective

#Agregamos una linea horizontal a la altura de M
abline(h=M,lty=2)
set.seed(3000)



#############################
#Distribucion MLE variando n#
#############################
iter=1000
res=NULL
N <- 7 #n
for(i in 1:iter){
#Definimos el vector de simulaciones
xsim <- numeric(0)
#Definimos la cantidad de simulaciones que queremos

#Esta es la rutina central de la simulación
while(length(xsim)<=(N-1)){
  x <- runif(1,min = -1,max=1)
  y <- runif(1,min=0,max=M)
  px <- p(x)
  if(y<=px){
    xsim <- c(xsim,x)
    color <- "green"
  }else{color <- "red"}
  
  points(x=x,y=y,pch=4,col=color)
  if(y<=px){
    points(x=x,y=0,pch=20)
  }
}
nlogL <- function(theta){ #loglikelihood negativo porque optim minimiza
  -sum(log(0.5*(1+theta*xsim)))
}
M0=optimize(f=nlogL,interval=c(-1,1))$minimum #estimador MLE
res=c(res,M0)
}

plot(density(res),main=paste('distribucion MLE_n=',N,', mean=',round(mean(res),3),', var=',round(var(res),2)))
abline(v=mean(res),lty=2)
legend(x="topleft",cex=0.7,horiz=F, c('mean'),col=c(1:2),lty=2,lwd=1)


    
#######################
#Información de Fisher#
#######################
rta <- c(2, 2, 1, 1, 1, 1, 0, 2, 1, 2,
         1, 0, 1, 2, 1, 0, 0, 2, 2, 1)
sum(dbinom(x=rta, size=5, prob=0.3, log=TRUE))
ll <- function(prob) sum(dbinom(x=rta, size=5, prob=prob, log=T))
plot(Vectorize(ll))
minusll <- function(x) -ll(x)
optimize(f=minusll, interval=c(0, 1))

#Informacion de Fisher
#Caso 1
w <- c(5)
ll1 <- function(lambda) sum(dpois(x=w, lambda=lambda, log=T))
ll1 <- Vectorize(ll1)
plot(ll1,xlim=c(0,15))
minusll1 <- function(x) -ll1(x)
optimize(f=minusll1, interval=c(0, 15))

# Caso 2
y <- c(5, 10, 6, 15)
ll2 <- function(lambda) sum(dpois(x=y, lambda=lambda, log=T))
ll2 <- Vectorize(ll2)
plot(ll2,xlim=c(0,15),add=T)
minusll2 <- function(x) -ll2(x)
optimize(f=minusll2, interval=c(0, 15))
# De la figura se observa claramente que cuando se tienen 4 observaciones 
# la curva es más puntiaguda y por lo tanto menor incertibumbre sobre el parámetro  
# λ a estimar.