# Lab 10 TSdD 2024 # Es. 1 # AfE - Set estimation - Normal model rm(list=ls()) mu.p=1.9 sig2.p=0.095 sd.p=sqrt(sig2.p) b=1 k=seq(0,4,.01) mis.S=2*sd.p*k prob.S=pnorm(k)-pnorm(-k) rho=b*mis.S-prob.S # (a) plot(k,rho,type="l") # (b) k.star=k[rho==min(rho)] k.star # verifica con risultato esatto k.star.esatto=sqrt(-2*log(sqrt(2*pi)*sd.p*b)) k.star.esatto # (c) L=mu.p-k.star*sd.p U=mu.p+k.star*sd.p C.star=c(L,U) C.star # (1.67, 2.12) # (d) mis.S.star=2*sd.p*k.star prob.S.star=pnorm(k.star)-pnorm(-k.star) rho.star=b*mis.S.star-prob.S.star rho.star # -0.08 # (e) # Ripeti con b=0.5 # Es. 2 # MC frequentista rm(list=ls()) M=10000 mu=3 sig2=2 sig=sqrt(sig2) x.mc=rnorm(M,3,sqrt(sig)) E.x=mean(x.mc) V.x=var(x.mc) y.mc=exp(x.mc) E.y=mean(y.mc) V.y=var(y.mc) E.x V.x E.y V.y # Es. 3 # MC Frequentist point estimation 1 # Unbiased estimator (Uniform model) rm(list=ls()) M=10000 n=10 theta.vero=0.4 x.matrix=runif(M*n,min=0,max=theta.vero) x.matrix=matrix(x.matrix,M,n) x.vett=apply(x.matrix,1,max) cc=(n+1)/n d.mv=cc*x.vett E.d.mv=mean(d.mv) E.d.mv var.d.mv=var(d.mv) var.d.mv var.vera=theta.vero^2/(n*(n+2)) var.vera # Es. 4. # MC frequentist point estimation 2 # Bias of a frequentist estimator (Beta model) rm(list=ls()) M=10000 n=10 theta.vero=2 x.matrix=rbeta(M*n,theta.vero,1) x.matrix=matrix(x.matrix,M,n) x.matrix=log(x.matrix) x.vett=apply(x.matrix,1,sum) x.vett=-n/x.vett E.d.mv=mean(x.vett) E.d.mv # Es. 5 rm(list=ls()) M=10000 n=3 R.fun=function(theta){ x.matr=runif(M*n,0,theta) x.matr=matrix(x.matr,M,n) d1=apply(x.matr,1,max) d2=((n+1)/n)*d1 d3=((n+2)/(n+1))*d1 B1=mean(d1)-theta V1=var(d1) R1=V1+B1^2 B2=mean(d2)-theta V2=var(d2) R2=V2+B2^2 B3=mean(d3)-theta V3=var(d3) R3=V3+B3^2 return(c(R1,R2,R3)) } R.fun(2) th.val=seq(0,5,0.1) th.val=matrix(th.val,length(th.val),1) R.val=apply(th.val,1,R.fun) R1=R.val[1,] R2=R.val[2,] R3=R.val[3,] plot(th.val,R1,type="l",xlab=expression(theta),lwd=2,ylim=c(0,2)) lines(th.val,R2,lty=2,lwd=2) lines(th.val,R3,lty=4,lwd=2)