# Lab 9 IS 2025 # Es. 1 # normalità asintotica S^2_n (UMVUE) modello N(theta1,theta2) rm(list = ls()) M=10000 n=10 theta1=1 theta2=4 # distribuzione empirica (MC) x.matrix=rnorm(M*n,theta1,sqrt(theta2)) x.matrix=matrix(x.matrix,M,n) S2.mc=apply(x.matrix,1,var) E.S2.mc=mean(S2.mc) V.S2.mc=var(S2.mc) E.S2.mc V.S2.mc # valore atteso e varianza esatti E.S2.esatta=theta2 V.S2.esatta=(2*theta2^2)/(n-1) E.S2.esatta V.S2.esatta # confronto istogramma (distribuzione empirica) / # distribuzione esatta di S^2_n / approssimazione normale a=(n-1)/2 b=(n-1)/(2*theta2) hist(S2.mc,freq=FALSE,ylim=c(0,0.25)) curve(dgamma(x,a,rate=b), add=TRUE, lty=2) # approssimazione normale m=theta2 v=2*theta2^2/n # cr(theta2) curve(dnorm(x,m,sqrt(v)), add=TRUE, lty=1) # ripeti per n=50 (ylim=c(0,0.6)) # Es. 2 # distribuzione SMV=1/X.med modello EN(theta) rm(list = ls()) M=10000 n=10 theta=2 # distribuzione empirica (MC) x.matrix=rgamma(n*M,1,rate=theta) x.matrix=matrix(x.matrix,M,n) medie=apply(x.matrix,1,mean) smv=1/medie E.smv.mc=mean(smv) V.smv.mc=var(smv) E.smv.mc V.smv.mc # E.smv.esatto=n*theta/(n-1) E.smv.esatto # distribuzione empirica (istogramma) - distribuzione esatta - approssimazione normale hist(smv,freq=F,ylim=c(0,0.8)) # usando funzione gamma inversa install.packages("invgamma") library(invgamma) help(dinvgamma) # distribuzione esatta curve(dinvgamma(x,n,rate=n*theta),lty=2,add=TRUE) # appross. normale m=theta v=theta^2/n curve(dnorm(x,m,sqrt(v)),lty=1,add=T) # ripeti per n=50 (ylim=c(0,1.6)) # Es. 3 # distribuzione di g(theta)=e^{-theta} modello Pois(theta) rm(list = ls()) M=10000 n=10 theta=2 # distribuzione empirica (MC) x.matrix=rpois(n*M,theta) x.matrix=matrix(x.matrix,M,n) smv=apply(x.matrix,1,mean) g=exp(-smv) E.g.mc=mean(g) V.g.mc=var(g) E.g.mc V.g.mc # appross. normali di E.g e V.g E.tilde.g=exp(-theta) V.tilde.g=exp(-2*theta)*theta/n E.tilde.g V.tilde.g # confronto distribuzione empirica e approx normale hist(g,freq=FALSE) curve(dnorm(x,E.tilde.g,sqrt(V.tilde.g)), add=TRUE) # per n=10 le appross. di MC ancora abbastanza imprecise # ripeti per n=50 # Es. 4 # IC normale N1, N2 e N3 (numerico) x.oss=c(63,65,94,37,83,95,70,96,47,29,52,38,47,79,66,25,48,80,52,49) n=length(x.oss) alpha=0.1 # (a) summary (x.oss) # (b) N1 -- IC per valore atteso - varianza nota varianza=400 n=length(x.oss) z = qnorm(1-alpha/2) m = mean(x.oss) L=m-z*sqrt(varianza/n) U=m+z*sqrt(varianza/n) c(L,U, U-L) # Per ottenere quanto chiesto: # confronta con install.packages("TeachingDemos") library(TeachingDemos) ? z.test z.test(x.oss,stdev=sqrt(400),conf.level=0.90) z.test(x.oss,stdev=sqrt(400),conf.level=0.90)[4] # (c) N2 -- IC per valore atteso - varianza incognita qt(1-alpha/2,n-1) L1=mean(x.oss)-qt(1-alpha/2,n-1)*sd(x.oss)/sqrt(n) U1=mean(x.oss)+qt(1-alpha/2,n-1)*sd(x.oss)/sqrt(n) IC1=c(L1,U1) IC1 # confronta con t.test(x.oss,conf.level=0.90)[4] # (d) N4 -- IC per varianza - valore atteso noto qchisq(1-alpha/2,n-1) qchisq(alpha/2,n-1) L2=(n-1)*var(x.oss)/qchisq(1-alpha/2,n-1) U2=(n-1)*var(x.oss)/qchisq(alpha/2,n-1) IC2=c(L2,U2) IC2 var(x.oss) # Esercizio 5 # Simulazione probabilità copertura normale N2 M=10000 n=10 alpha=0.2 th.vero=-3 sigma2=1 # (a) x.matr=matrix(rnorm(M*n,th.vero,sqrt(sigma2)),M,n) # (b) medie=apply(x.matr,1,mean) S.n=apply(x.matr,1,sd) t.al=qt(1-alpha/2,n-1) C.inf=medie-t.al*S.n/sqrt(n) C.sup=medie+t.al*S.n/sqrt(n) C.inf[1:10] C.sup[1:10] # (c) mean(C.sup<=2.5) mean(C.inf>2.5) # (d) copertura=mean(C.infth.vero) copertura # ES 6 # Confronto ET e HD modello EN n=10 alpha=0.05 # ET interval q1.et=qgamma(alpha/2,n,1) q2.et=qchisq(1-alpha/2,n,1) q1.et q2.et L.et=q1.et/T U.et=q2.et/T C.et=c(L.et,U.et) Lungh.et=U.et-L.et C.et Lungh.et Lungh.hd