prova di esame

1.Date U1,U2 … iid U(0,1), si descriva il metodo dell’inversione della funzione di ripartizione per generare delle determinazioni casuali da una variabile casuale (v.c.) Y continua con f(y) = 3y2 per 0<y<1

cdf <- function(y) {y^3}
inverse.F <- function(u) {u^(1/3)}
R <- 1000
# campiono R valori da U=F(Y) ~ U(0,1)
u <- runif(R)
# calcolo y = F-1(u) 
y <- inverse.F(u=u)

2.Si descriva il metodo dell’inversione della funzione di ripartizione per generare delle determinazioni casuali da una variabile casuale (v.c.) Y discreta che assume valore 1 con probabilità 0.6 e valore 0 con probabilità 0.4.

p = 0.6
supporto = 0:1
Pi = pbinom(q=0:1, size=1, prob=p)

# campiono da U(0,1) e fisso y = min(k: Pk>=u)
u = runif(1000)
y = c()
for (i in 1:1000) {
  y[i] = supporto[min(which(Pi>=u[i]))]
}

3. Si descriva il metodo di accettazione-rifiuto per generare delle determinazioni dalla v.c. Y continua con funzione di densità: f(y) = 3y2 con (0 < y < 1), scegliendo come gx(.) la funzione di densità di una Uniforme in (0,1)
   

# k = max(f/g) = maxf
f <- function(y) {3*y^2}
k <- optimize(f, lower=0, upper=1, maximum=T)$objective
R <- 1000
u <- runif(R)
x <- runif(R)
# test: u <= f/(k*g) = f/k
test <- (u<=(f(y=x)/k))
# campione di valori accettati
y <- x[test==T]

4. Si generi un campione di ampiezza n = 30 determinazioni pseudo-casuali da una v.c. Normale con media µ = 2 e varianza σ2 = 5, avendo cura di fissare il seme dell’estrazione tramite il comando set.seed(345).
   
5. Trovare la stima di massima verosimiglianza (SMV) per il vettore dei parametri (µ,σ) attraverso il comando optim(), usando come punti iniziali la mediana campionaria e la radice quadrata del momento secondo campionario.
   
6. Individuare un intervallo di confidenza alla Wald per entrambi i parametri con livello di confidenza 0.98.
   
7. Trovare la SMV per lo stesso modello utilizzando la ri-parametrizzazione (µ,ζ) = (µ,1/σ).

n <- 30
media <- 2
varianza <- 5
set.seed(345)

y <- rnorm(n, mean=media, sd=sqrt(varianza))

mediana <- median(y)
sqrt.x2 <- sqrt(mean(y^2))
start <- c(mediana, sqrt.x2)

nlogL <- function(param, data) {
-sum(dnorm(data, mean=param[1], sd=sqrt(param[2]), log=T))
}

#Vectorize(nlogL, vectorize.args='param')

ML <- optim(start, fn=nlogL, data=y, hessian=T)
stima <- ML$par
stima

## [1] 1.852404 5.721796

var.stima <- solve(ML$hessian)
var.stima

##              [,1]         [,2]
## [1,] 1.907265e-01 9.134125e-05
## [2,] 9.134125e-05 2.182389e+00

conf.level = 0.98
alpha = 1 - conf.level
ic.wald.media <- c(stima[1] - qnorm(1-alpha/2)*sqrt(var.stima[1,1]), stima[1] + qnorm(1-alpha/2)*sqrt(var.stima[1,1]))

ic.wald.var <- c(stima[2] - qnorm(1-alpha/2)*sqrt(var.stima[2,2]),
                 stima[2] + qnorm(1-alpha/2)*sqrt(var.stima[2,2]))

# sfrutto proprietà di invarianza dello SMV 
# se ho SMV per * e ne voglio uno per f(*), applico direttamente la trasformazione f
ic.wald.eta <- c(1/stima[2] - qnorm(1-alpha/2)*sqrt(1/var.stima[2,2]),
                 1/stima[2] + qnorm(1-alpha/2)*sqrt(1/var.stima[2,2]))