# Lab 4 IS 2025 -- 20 marzo 2025


# Es. 1 (fdv modello N1)
# Funzione per fdv modello N1
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. 2 (insiemi Lq per esempio N1)
# codice (non in forma di funzione) per insiemii Lq modello N1
dati1=c(-2,3,4,2,6)
dati2=c(dati1,dati1,dati1)
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
# Funzione per insiemi di verosimiglianza (esatti) modello N1
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. 3 (insiemi Lq per esempio N1: numericamente)
# codice (non in forma di funzione) per insiemii Lq
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



# Es. 4 (modello Bernoulliano: approssimazione normale di L)
# funzione per fdv relativa esatta
L.bern.rel=function(theta,xn){
  sn=sum(xn)
  n=length(xn)
  smv=sn/n
 N=(theta^sn)*((1-theta)^(n-sn))
 D=(smv^sn)*((1-smv)^(n-sn))
 return(N/D)
 # NB: N/D si può ottenere anche come segue
 # dbeta(theta,shape1=sn+1,shape2=n-sn+1)/dbeta(smv,shape1=sn+1,shape2=n-sn+1)
  }
# funzione per fdv relativa, approssimazione normale 
L.bern.approx=function(theta,xn){
  smv=mean(xn)
  n=length(xn)
  I.oss=n/(smv*(1-smv))
  A=I.oss/2
  B=(theta-smv)^2
  exp(-A*B)
}
dati1=c(rep(0,5),rep(1,10))
dati2=c(rep(0,15),rep(1,30))
curve(L.bern.rel(x,xn=dati1),from=0,to=1,xlab=expression(theta),ylab=expression(L(theta)))
curve(L.bern.approx(x,xn=dati2),add=TRUE,lty=2)



# Es. 5 (modello Bernoulliano: insiemi Lq approssimati)
# funzione per intervallo esatto determinato numericamente
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)
abline(h=0.5)
abline(v=0.517)
abline(v=0.796)
Lq.bern.num(xn=dati2,q=0.5)
L=Lq.bern.num(xn=dati2,q=0.5)[1]
U=Lq.bern.num(xn=dati2,q=0.5)[2]
L
U
abline(v=L,lty=2)
abline(v=U,lty=2)
# funzione per intervallo approssimato
Lq.bern.appr=function(xn,q){
q=0.5
n=length(xn)
smv=mean(xn)
I.oss=n/(smv*(1-smv))
k=sqrt(-2*log(q))
Ln.appr=smv-k*sqrt(1/I.oss)
Un.appr=smv+k*sqrt(1/I.oss)
Lungh.appr=Un.appr-Ln.appr
return(c(Ln.appr,Un.appr,Lungh.appr))
}
Lq.bern.num(xn=dati2,q=0.5)
Lq.bern.appr(xn=dati2,q=0.5)
# altri valori di q
Lq.bern.num(xn=dati1,q=0.147)
Lq.bern.appr(xn=dati1,q=0.147)
# con pochissime osservazioni
dati0=c(0,1,1)
Lq.bern.num(xn=dati0,q=0.5)
Lq.bern.appr(xn=dati0,q=0.5)