# 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)