# Lab. 5 IS 2025

# Esercizio 1

# fdv esatta (non relativa)
L.ex=function(theta,x,mu){
  n=length(x)
  S02=mean((x-mu)^2)
  (theta^(-n/2))*exp(-(n*S02)/(2*theta))
}

# fdv relativa esatta
L.rel=function(theta,x,mu){
  S02=mean((x-mu)^2)
  L.ex(theta,x,mu)/L.ex(S02,x,mu)
}

# Approssimazione normale fdv
L.appr=function(theta,x,mu){
  n=length(x)
  smv=mean((x-mu)^2)
  I.n=n/(2*(smv^2))
  exp(-0.5*((theta-smv)^2)*I.n)
}

# (b)
dati=c(2.5,0.8,2.7,2.4,4.1,6.1,1,3,2.9,-2)
mu.0=2
curve(L.rel(x,dati,mu.0), from=0, to=20, xlab=expression(theta), ylab="verosimiglianza relativa")

# (c)
curve(L.appr(x,dati,mu.0),add=TRUE,lty=2)

# (d)
dati2=c(dati,3,5,6,1,1)
curve(L.rel(x,dati2,mu.0), from=0, to=20, xlab=expression(theta), ylab="verosimiglianza relativa")
curve(L.appr(x,dati2,mu.0),add=TRUE,lty=2)

# (e) Lq appross
Lq.appr=function(x,mu,q){
  n=length(x)
  smv=mean((x-mu)^2)
  I.n=n/(2*(smv^2))
  kq=sqrt(-2*log(q))
  L.n=smv-kq*sqrt(1/I.n)
  U.n=smv+kq*sqrt(1/I.n)
  return(c(smv,L.n,U.n))
  }

Lq.appr(dati,mu=2,q=0.8)



# Esercizio 2

# (a)
L.fun=function(theta,x){
  n=length(x)
  0*(theta<max(x)) + (theta^(-n))*(theta>=max(x))
}

# (b) 
L.fun.rel=function(theta,x){
  L.fun(theta,x)/L.fun(max(x),x)
}

dati=c(4.6, 1.0, 1.0, 0.9, 3.1, 1.7, 4.7, 0.6, 4.0, 1.3)
# SMV
max(dati)

curve(L.fun.rel(x,dati),from=0,to=10,ylab="fdv relativa ", xlab=expression(theta),ylim=c(0,1))
title("Funzione di verosimiglianza relativa -- modello uniforme")

# Con “abline()” disegniamo la retta y=0.147
abline(h=0.147)

# (c) Insieme L.q
# Per trovare l'insieme L.q, definiamo una griglia di valori
# (ad esempio tra 2 ed 8 di passo 0.0001)
th.val=seq(2,8,.0001)

# Troviamo l'estremo inferione dell'intervallo (th.L) selezionando
# il più piccolo valore del vettore th.val in corrispondenza del quale
# la funzione L.fun.rel() risulta superiore a 0.147 ---> OVVIAMENTE viene x_(n)=max(x_1,...,x_n)
th.L=min(th.val[L.fun.rel(th.val,dati)>=0.147])
th.L

# Analogamente, l'estremo superiore dell'intervallo (th.U)
# lo troviamo selezionando il più grande valore del vettore in corrispondenza del quale L.fun.rel()>0.147
th.U=max(th.val[L.fun.rel(th.val,dati)>=0.147])
th.U

# L'intervallo cercato è quindi ottenuto concatenando
# (con la funzione “c()”) i due valori determinati, th.L e th.U
c(th.L,th.U)

# (4) La lunghezza dell'intervallo si ottiene calcolando
# la differenza tra gli estremi dell'intervallo
Lungh=th.U-th.L
Lungh

# (5) Per ottenere il nuovo vettore dei dati basta concatenare “dati1”
# alle nuove osservazioni campionarie
# n=15
dati2=c(dati,3.7,4.2,2.1,4.3,0.9)

# Usiamo quindi “curve” per aggiungere il nuovo grafico, ricordando di
# specificare “add=T” per sovrapporre il nuovo al vecchio grafico
curve(L.fun.rel(x,dati2),lty=2,add=T)

# (6) Per trovare il nuovo insieme L.q ripetiamo quanto fatto
# al punto (3)
# Insieme L.q
th.val=seq(2,8,.0001)
th.L=min(th.val[L.fun.rel(th.val,dati2)>=0.147])
th.L
th.U=max(th.val[L.fun.rel(th.val,dati2)>=0.147])
th.U
c(th.L,th.U)
Lungh=th.U-th.L
Lungh

# (7) da soli per casa



# Es 3
thetavero=3
sig2=1
sig=sqrt(1)
n=1
curve(dnorm(x,mean=thetavero,sd=sig),lwd=2,from=-1,to=7,
      xlab="valori dello stimatore Xmed",
      ylim=c(0,5),ylab="funzione di densità dello stimatore")
abline(v=thetavero)
n=10
curve(dnorm(x,mean=thetavero,sd=sig/sqrt(n)),add=TRUE,lty=2)
n=50
curve(dnorm(x,mean=thetavero,sd=sig/sqrt(n)),add=TRUE,lty=3)
n=100
curve(dnorm(x,mean=thetavero,sd=sig/sqrt(n)),add=TRUE,lty=4)
# n=1000
# curve(dnorm(x,mean=thetavero,sd=sig/sqrt(n)),add=TRUE,lty=5)



# Es 4
thetavero=2
n=5
curve(dgamma(x,n/2,rate=n/(2*thetavero)),lwd=2,from=0,to=10,xlab="S02",
      ylim=c(0,2),ylab="funzione di densità dello stimatore")
abline(v=thetavero)
n=10
curve(dgamma(x,n/2,rate=n/(2*thetavero)),add=TRUE,lty=2)
n=50
curve(dgamma(x,n/2,rate=n/(2*thetavero)),add=TRUE,lty=3)
n=100
curve(dgamma(x,n/2,rate=n/(2*thetavero)),add=TRUE,lty=4)
n=1000
curve(dgamma(x,n/2,rate=n/(2*thetavero)),add=TRUE,lty=5)