# LIS 24-25
# Laboratorio 3 - 13 marzo 2025
# Funzione  di verosimiglianza (fdv)


# Esercizio 1 (modello Bernoulliano: fdv)
rm(list = ls())
L.bern=function(x,theta){
  s=sum(x)
  n=length(x)
  (theta^s)*((1-theta)^(n-s))
}
dati=c(1,1,0,0,1)
curve(L.bern(dati,x),from=0,to=1, xlab=expression(theta),ylab="verosimiglianza", lwd=2)
# oppure
# th.val=seq(0,1,0.01)
# plot(th.val,L.bern(dati,th.val),type="l",lwd=2,ylab="verosimiglianza",xlab=expression(theta))
# Punto di massimo di della fdv
th.val=seq(0,1,0.01)
th.max=th.val[which.max(L.bern(dati,th.val))]
th.max
# controllo
mean(dati)



# Es. 2 (modello Bernoulliano: fdv relativa)
dati1=c(rep(0,5),rep(1,10))
dati2=c(rep(0,15),rep(1,30))
L.bern.rel=function(theta,xn){
  sn=sum(xn)
  n=length(xn)
  smv=sn/n
  dbeta(theta,shape1=sn+1,shape2=n-sn+1)/dbeta(smv,shape1=sn+1,shape2=n-sn+1)
}
curve(L.bern.rel(x,xn=dati1),from=0,to=1,xlab=expression(theta),ylab=expression(L(theta)))
curve(L.bern.rel(x,xn=dati2),add=TRUE,lty=2)



# Es. 3 (modello Bernoulliano: Lq numerico)
Lq.bern.num=function(xn,q){
  th.value=seq(0,1,.001)
  L.value=L.bern.rel(th.value,xn)
  # estremo inferiore intervallo: L--> lower
  Ln.num=min(th.value[L.value>q])
  # estremo superiore intervallo: U--> upper
  Un.num=max(th.value[L.value>q])
  Lungh=Un.num-Ln.num
  return(c(Ln.num,Un.num,Lungh))
}
Lq.bern.num(xn=dati1,q=0.5)
# 0.517, 0.796, 0.27
abline(h=0.5)
abline(v=0.517)
abline(v=0.796)
# con n più grande
Lq.bern.num(xn=dati2,q=0.5)
abline(v=0.582,lty=2)
abline(v=0.745,lty=2)



# Es. 4 (fdv modello N1)
rm(list = ls())
L.N1=function(theta,xn,sig2){
  xmed=mean(xn)
  n=length(xn)
  A=n/(2*sig2)
  B=(theta-xmed)^2
  exp(-A*B)
}
dati1=c(-2,3,4,2,6)
mean(dati1)
length(dati1)
curve(L.N1(x,xn=dati1,sig2=1), from=0, to=5, ylab=expression(L(theta)), xlab=expression(theta))
dati2=c(dati1,dati1,dati1)
mean(dati2)
length(dati2)
curve(L.N1(x,xn=dati2,sig2=1),add=TRUE,lty=2)


# Es. 5 (insiemi Lq per esempio N1)
xn=dati1
n=length(xn)
xmed=mean(xn)
sig2=1
sig=sqrt(sig2)
q=0.8
k=sqrt(-2*log(q))
Ln=xmed-k*sig/sqrt(n)
Un=xmed+k*sig/sqrt(n)
c(Ln,Un)
Lungh=Un-Ln
Lungh
# come funzione
Lq.N1=function(xn,sig2,q){
  n=length(xn)
  xmed=mean(xn)
  sig=sqrt(sig2)
  k=sqrt(-2*log(q))
  Ln=xmed-k*sig/sqrt(n)
  Un=xmed+k*sig/sqrt(n)
  Lq=c(Ln,Un)
  Lungh=Un-Ln
  return(c(Lq,Lungh))
}
Lq.N1(xn=dati1,sig2=1,q=0.8)
Lq.N1(xn=dati2,sig2=1,q=0.8)


# Es. 6 (insiemi Lq per esempio N1: numericamente)
th.value=seq(1,5,.01)
L.value=L.N1(th.value,xn=dati1,sig2=1)
q=0.8
# estremo inferiore intervallo: L--> lower
Ln.num=min(th.value[L.value>q])
Ln.num
abline(v=Ln.num)
# estremo superiore intervallo: U--> upper
Un.num=max(th.value[L.value>q])
Un.num
abline(v=Un.num)
abline(h=q)
c(Ln.num,Un.num)
# Ripeti con dati2