# Lab 7 - 10 aprile 2025


# Es 1 (vedi esempio 4.28 dispense)
# Confronto MSE smv e smo (modello uniforme)
MSE.smv=function(theta,n){
  theta^2/((n+1)*(n+2))
}
MSE.smo=function(theta,n){
  theta^2/(3*n)
}
curve(MSE.smv(x,20),from=0,to=10,
      ylim=c(0,1),lwd=2,xlab=expression(theta),ylab="MSE")
curve(MSE.smo(x,20),lwd=2,lty=2, add=TRUE)
# Efficienza
MSE.smo(theta=6,n=20)/MSE.smv(theta=6,n=20)
# Es 2 (vedi esempio 4.69 dispense)
# Confronto MSE varianza campionaria (vc) e 
# varianza campionaria corretta (vcc) (modello normale)
MSE.vcc=function(theta,n){
  2*theta^2/(n-1)
}
MSE.vc=function(theta,n){
(2*n-1)*theta^2/(n^2)
}
curve(MSE.vcc(x,5),from=0,to=5,
      ylim=c(0,15),lwd=2,xlab=expression(theta),ylab="MSE")
curve(MSE.vc(x,5),lwd=2,lty=2, add=TRUE)



# Es 3 (vedi esempio 4.30 dispense)
# Modello bernoulliano: 
# (a) consistenza smv 
# (b) confronto con altro stimatore (con MSE)
# Consistenza
MSE.Xmed=function(theta,n) { theta*(1-theta)/n }
curve(MSE.Xmed(x,10),from=0,to=1,
      ylim=c(0,0.03),lwd=2,xlab=expression(theta),ylab="MSE")
curve(MSE.Xmed(x,30),add=T,lty=2,lwd=2)
curve(MSE.Xmed(x,50),add=T,lty=3,lwd=2)
# per qualsiasi theta fissato, ad esempio theta=0.4
MSE.Xmed(theta=0.4,n=1:100)
plot(1:100,MSE.Xmed(theta=0.4,n=1:100),type="l")
# Confronto con MSE del secondo stimatore d2
curve(MSE.Xmed(x,10),from=0,to=1,
      ylim=c(0,0.03),lwd=2,xlab=expression(theta),ylab="MSE")
n=10
abline(h=n/(4*(n+sqrt(n))^2),lwd=2)
a=0.5-0.5*sqrt(1-(n^2/(n+sqrt(n))^2))
abline(v=a)           
b=0.5+0.5*sqrt(1-(n^2/(n+sqrt(n))^2))
abline(v=b)           
c(a,b)



# Es 4 (chi quadrato)
# Chi quadrato
dchisq(x,df) # funzione di densità
pchisq(x,df) # f.ne ripartizione
qchisq(alpha,df) # quantili
rchisq(M,df=) # generatore n. pseudocasuali
# (a) densità
# dchisq(x,df)
dchisq(0.5,3)
curve(dchisq(x,df=1),from=0,to=20,lwd=2,ylab="densità chi quadrato")
curve(dchisq(x,df=3),from=0,to=20,lwd=2,add=TRUE,lty=2)
curve(dchisq(x,df=5),from=0,to=20,lwd=2,add=TRUE,lty=3)
curve(dchisq(x,df=10),from=0,to=20,lwd=2,add=TRUE,lty=4)
# (b)
gdl=10
curve(dchisq(x,df=gdl),from=0,to=4*gdl,lwd=2,ylab="densità chi quadrato")
curve(dnorm(x,mean=gdl,sd=sqrt(2*gdl)),lwd=2,add=TRUE,lty=2)
# (c)
gdl=30
curve(dchisq(x,df=gdl),from=0,to=60,lwd=2,ylab="densità chi quadrato")
curve(dnorm(x,mean=gdl,sd=sqrt(2*gdl)),lwd=2,add=TRUE,lty=2)



# Es 5 -- Binom(n,theta)
n=400
theta=0.1
# (a)
sum(dbinom(35:n,size=n,prob=theta))
# (b)
sum(dbinom(0:20,size=n,prob=theta))
# (c)
sum(dbinom(0:1,size=n,prob=theta))
# (d)
sum(dbinom(40:50,size=n,prob=theta))



# Es 6 -- Binom (n,theta)
n=100
theta=0.035
sum(dbinom(3:n,size=n,prob=theta))



# Es 7 -- Normale (media campionaria)
n=9
mu=12.2
sig=3.6
1-pnorm(14.4,mu,sig/sqrt(n))



# Es 8 -- Normale (somma di somme campionarie indipendenti)
nA=20
nB=30
muA=25
muB=40
sig2A=64
sig2B=81
# Y=valore portafoglio 
# Y~N(muY,sig2Y)
muY=nA*muA+nB*muB
sig2Y=sig2A*nA+sig2B*nB
# (a)
1-pnorm(1800,muY,sqrt(sig2Y))f
# 
pnorm(1650,muY,sqrt(sig2Y))