# LIS 23-24 CASI DI STUDIO

# Normal model (gasoline consumption, Newbold et al (2010) p. 325)
# consumo medio settimanale per miglia di gasolio per 24 furgoni
# la varianza la fisso io: sig2=2.8 (circa uguale a S2_n )

# (a)
gasoline=scan()
15.5 
21.0 
18.5 
19.3 
19.7 
16.9 
20.2 
14.5
16.5 
19.2 
18.7 
18.2 
18.0 
17.5 
18.5 
20.5
18.6 
19.1 
19.8 
18.0 
19.8 
18.2 
20.3 
21.8

sigma2=2.8   # considerata nota, dato del problema


# (b)
summary(gasoline)
# (c)
hist(gasoline,freq=FALSE,ylim=c(0,0.4))
# (d)
lines(density(gasoline))
# (e)
mean(gasoline)
# (f)
curve(dnorm(x,mean(gasoline,sqrt(sigma2))),add=TRUE,lty=2)
# vedi anche
shapiro.test(gasoline)
qqnorm(gasoline)
qqline(gasoline)
boxplot(gasoline)

# (g) 
# fdv relativa [funzione generale]
Lik.norm.rel=function(theta,y,sig2){
  y.med=mean(y)
  n=length(y)
  sig2=2.8
  A=n*((theta-y.med)^2)/(2*sig2)
  exp(-A)
}
curve(Lik.norm.rel(x,gasoline,2.8),from = 17, to=20,xlab=expression(theta),ylab="fvd relativa")
smv=mean(gasoline)
abline(v=smv)


# (h) 
# insieme L.q
L.q.fun=function(q,y,sig2){
k.q=sqrt(-2*log(q))
smv=mean(y)
n=length(y)
sig=sqrt(sig2)
th.L=smv-k.q*sig/sqrt(n)
th.U=smv+k.q*sig/sqrt(n)
L.q=c(th.L,th.U)
return(L.q)
}
L.q.fun(0.75,gasoline,2.8)


# (i) 
# verifica grafica
th.L.gasoline=L.q.fun(0.75,gasoline,2.8)[1]
th.U.gasoline=L.q.fun(0.75,gasoline,2.8)[2]
abline(h=0.75)
abline(v=th.L.gasoline)
abline(v=th.U.gasoline)


# (j)
# Altro modo di procedere: scriviamo le due funzioni 
# in funzione di (n,smv,sig2)
# Per la fdv relativa
Lik.norm.rel.SMV=function(theta,n,smv,sig2){
  A=n*((theta-smv)^2)/(2*sig2)
  exp(-A)
}
curve(Lik.norm.rel.SMV(x,n=length(gasoline),smv=mean(gasoline),sig2=2.8),
      from = 17, to=20,xlab=expression(theta),ylab="fvd relativa")
# vediamo se consideriamo un campione con stessa media di gasoline ma n=100
curve(Lik.norm.rel.SMV(x,n=100,smv=mean(gasoline),sig2=2.8), lty=2,add=TRUE)


# (k)
# Per la funzione per L.q
L.q.fun.SMV=function(q,n,smv,sig2){
  sig=sqrt(sig2)
  k.q=sqrt(-2*log(q))
  th.L=smv-k.q*sig/sqrt(n)
  th.U=smv+k.q*sig/sqrt(n)
  L.q=c(th.L,th.U)
  return(L.q)
}
# per n=24, smv=mean(gasoline)
L.q.fun.SMV(q=0.75,n=length(gasoline),smv=mean(gasoline),sig2=2.8)
# per n=100, smv=mean(gasoline)
L.q.fun.SMV(q=0.75,n=100,smv=mean(gasoline),sig2=2.8)