# Lab 8 IS 2025 -- 24 aprile 2022


# Esercizio 1
M=10000
alpha=2
beta= 3
x.simul=rgamma(M,alpha,rate=beta)
hist(x.simul,prob=T)
curve(dgamma(x,alpha,rate=beta),add=T)
# (a) valore atteso simulato
EX=mean(x.simul)
EX
# valore atteso esatto
alpha/beta
# (b) varianza simulata
VX=var(x.simul)
VX
# varianza esatta
alpha/beta^2
# (c) probabilità simulata
mean(x.simul>0.5&x.simul<2)
# probabilità esatta
pgamma(2,alpha,rate=beta)-pgamma(0.5,alpha,rate=beta)
# (d) mediana simulata
quantile(x.simul,0.7)
# mediana esatta
qgamma(0.7,alpha,rate=beta)
# (e) quantile simulata
quantile(x.simul,0.7)
# quantiile esatto
qgamma(0.7,alpha,rate=beta)



# Esercizio 2
install.packages("invgamma")
library(invgamma)
M=50000
alpha=4
beta= 3
x.simul=rgamma(M,alpha,rate=beta)
y.simul=1/x.simul
hist(y.simul,prob=T)
# (a) Valore atteso
EY=mean(y.simul)
EY
# valore atteso esatto: 
beta/(alpha-1)
# (b) Varianza
VY=var(y.simul)
VY
# varianza esatta: 
(beta^2)/(((alpha-1)^2)*(alpha-2))
# (c) Probab intervallo
mean(y.simul>0.5&y.simul<2)
# esatto
pinvgamma(2,alpha,rate=beta)-pinvgamma(0.5,alpha,rate=beta)
# (d) - (e) Mediana e quantile
quantile(y.simul,0.5)
quantile(y.simul,0.7)
# esatti
qinvgamma(0.5,alpha,rate=beta)
qinvgamma(0.7,alpha,rate=beta)



# Esercizio 3
M=10000
a=1
b=0.5
x.simul=rbeta(M,shape1=a,shape2=b)
hist(x.simul,prob=T)
curve(dbeta(x,a,b),add=T)
y.simul=-log(x.simul)
hist(y.simul,prob=T)
# (a) E[X] e E[Y]
E.X=mean(x.simul)
E.Y=mean(y.simul)
E.X
E.Y
# E[X]=a/(a+b)
a/(a+b)
# (b) V[X] e V[Y]
V.X=var(x.simul)
V.Y=var(y.simul)
V.X
V.Y
# Varianza esatta beta... vedi Wiki
# (c) P[X>0.2] e P[Y>0.2]
# Per X abbiamo
1-pbeta(0.2,a,b)
# con dati simulati
mean(x.simul>0.2)
# oppure
length(x.simul[x.simul>0.2])/M
# Per Y abbiamo
mean(y.simul>0.2)
# (d) P[0.1<X<0.3] e P[0.1<Y<0.3]
# Per X
pbeta(0.3,a,b)-pbeta(0.1,a,b)
# oppure
mean(x.simul<0.3&x.simul>0.1)
# Per Y
mean(y.simul<0.3&y.simul>0.1)
# (e) Mediane
# Per X
qbeta(0.5,a,b)
# oppure
median(x.simul)
# Per Y
median(y.simul)
# verifico
mean(y.simul<median(y.simul))
# (f) Quantili 0.9
# Per X
qbeta(0.9,a,b)
# oppure
quantile(x.simul,0.9)
# Per Y
quantile(y.simul,0.9)



# Esercizio 4
M=10000
n=10
mu=2
sig2=4
sig=sqrt(sig2)
# (a)
x.matrice=rnorm(n*M,mean=mu, sd=sig)
x.matrice=matrix(x.matrice,M,n)
# (b)
x.vett=apply(x.matrice,1,mean)
# (c)
mean(x.vett)
# (d)
var(x.vett)
# varianza esatta di X.med
4/n
# (e)
sig2=3
sig=sqrt(sig2)
x.matrice=rnorm(n*M,mean=1, sd=sig)     # e poi ripeti quanto fatto in (a)-(d)
# (f)
sig2=4
sig=sqrt(sig2)
hist(x.vett, ylim=c(0,0.8), freq=F)
curve(dnorm(x,2,sig/sqrt(n)),-1,5, add=T)



# Esercizio 5
M=10000
n=10
alpha=2
beta=3
x.matrice=rgamma(M*n,alpha,rate=beta)
x.matrice=matrix(x.matrice,M,n)
t.simul=apply(x.matrice,1,mean)
ET=mean(t.simul)
alpha/beta
VT=var(t.simul)
VT
(alpha/beta^2)/n



# Esercizio 6
M=10000 
n=10
sig2=4
sig=sqrt(sig2)
# (a)
x.matrice=rnorm(n*M,mean=2, sd=sig)
x.matrice=matrix(x.matrice,M,n)
# (b)
x.vett=apply(x.matrice,1,var)
# passo da S2.n (varianza campionaria corretta) a hat.sigma2.n (varianza campionaria non corretta)
x.vett=(n-1)*x.vett/n
# (c)
mean(x.vett)
# controllo
((n-1)/n)*sig2
# (d)
var(x.vett)
# controllo
2*(n-1)*sig2^2/(n^2)
# (e)
sig2=3
sig=sqrt(sig2)
x.matrice=rnorm(n*M,mean=1, sd=sig) # e poi ripeti
# (f)
var.fun=function(n){
  M=1000
  m=2
  sig=2
  x.matrice=matrix(rnorm(M*n,mean=m,sd=sig),M,n)
  x.vett=(n-1)*apply(x.matrice,1,var)/n
  return(var(x.vett))
}
# Per disegnare il grafico (ad esempio per n tra 2 e 100)
nn=seq(2,100,1)
nn=matrix(nn,length(nn),1)
var.val=apply(nn,1,var.fun)
plot(nn,var.val,type="l")
# controllo con varianza esatta dello stimatore
sigma=2
var.vera=2*(nn-1)/(nn^2)*(sigma^4)
lines(nn,var.vera,lty=2)
plot(nn,var.vera,type="l",lty=2)