Metodi Monte Carlo (MC)

Abbiamo una v.c. Y con funzione di probabilità g(y).
L’obiettivo è calcolare la media di una sua trasformazione: E[m(Y)]
Anziché calcolare l’integrale nel supporto di Y di m(y)g(y) in modo analitico, ricorriamo alla simulazione Monte Carlo: si generano R valori y da g() e si calcola la media empirica 1/R*sum(m(y)).

# esempio 1
# m(y) = (cos(50y)+sin(20y))^2
# y idd U(0,1)
require(ggplot2)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.6.3
R <- 1000
y <- runif(R)
m <- function(x) {(cos(50*x)+sin(20*x))^2}

# integrale analitico
integrale <- integrate(m, lower=0, upper=1)

# metodo MC
stima.MC <- mean(m(y))

# confronto
c(stima.MC, integrale$value)
## [1] 0.9199768 0.9652009
# stima varianza stimatore MC
var.MC <- (1/R)*(1/(R-1))*sum((m(y)-stima.MC)^2)

# IC
c(stima.MC - 2*sqrt(var.MC), stima.MC + 2*sqrt(var.MC))
## [1] 0.8543844 0.9855692
# rappresentazione grafica
plot(x=sort(y),y=m(sort(y)), type='l')

my <- m(y)
media.my <- c()
media.my[1] <- my[1]
for (i in 2:R) {
  media.my[i] <- mean(my[1:i])
}
upper.limit <- media.my + 2*sqrt(var.MC)
lower.limit <- media.my - 2*sqrt(var.MC)
df = data.frame(x=1:R, y=media.my, a=upper.limit, b=lower.limit)
ggplot() +
  geom_line(data=df, aes(x=x,y=y), color='blue') +
  geom_line(data=df, aes(x=x,y=a), color='red') +
  geom_line(data=df, aes(x=x,y=b), color='red')

Lo stesso principio può essere applicato per calcolare, anziché il valore atteso di una trasformazione m(Y), anche per calcolare la probabilità di eventi. Ad esempio Pr(Y<=t) si può simulare generando R valori da g() e calcolando la proporzione di valori simulati <=t.

# Y con pdf N(0,1)
# stima MC di P(Y<1)
R <- 5000
y <- rnorm(R,0,1)
stima.MC <- sum(y<1)/length(y)
prob.esatta <- integrate(dnorm,lower=-Inf,upper=1)$value
# confronto 
c(stima.MC, prob.esatta)
## [1] 0.8404000 0.8413448

Esercizio: simulando numeri casuali da U(0,1) calcola con metodo MC l’integrale definito da 1 a 3 di y/(1+y2)2

R <- 10000
u <- runif(R)
f <- function(y) {y/((1+y^2)^2)}
# voglio y in (1,3) anziché (0,1)
# se U è unif(0,1) -> Y = (b-a)*U + a è unif(a,b)
y <- 2*u + 1
hist(y, freq=F)

int.esatto <- integrate(f, lower=1, upper=3)$value
m <- function(y)  {2*f(y)}
stima.mc <- sum(m(y))/R
c(int.esatto,stima.mc)
## [1] 0.2000000 0.1986535
# IC 
var.mc <- sum((m(y)-stima.mc)^2)/(R*(R-1))
upper <- stima.mc + 2*sqrt(var.mc)
lower <- stima.mc - 2*sqrt(var.mc)
IC <- c(lower,upper); IC
## [1] 0.1961918 0.2011151
mc <- c()
mc[1] <- m(y[1])
for (i in 2:R)  {
  mc[i] <- mean(m(y[1:i]))
}
upper <- mc + 1.64*sqrt(var.mc)
lower <- mc - 1.64*sqrt(var.mc)
df <- data.frame(x=1:R, y=mc, a=upper, b=lower)
ggplot() +
  geom_line(data=df, aes(x=x, y=y), color='blue') +
  geom_line(data=df, aes(x=x, y=a), color='red') +
  geom_line(data=df, aes(x=x, y=b), color='red')

Proprietà dello stimatore MC

1 Lo stimatore MC per la media della trasformata m(Y) è non distorto e, per la legge dei grandi numeri, consistente.
2 La varianza dello stimatore è pari a var(m(Y))/R e si può stimare utilizzando gli R valori y simulati.
3 Per il TLC si possono costruire IC di livello alpha per la media di m(Y) come: stima.MC +- quantilesqrt(stima.var.MC)
4 Si può fissare l’accuratezza dello stimatore (ad es. al terzo decimale) risolvendo per R la disequazione: quantilesqrt(stima.var.MC) < 10^(-3)

Esercizio: data Z v.c. N(0,1), calcola P(Z>20) con metodo MC 

R <- 10^8
set.seed(123)
# f densità N(0,1) da integrare da 20 a +Inf
# cambio variabile z = 1/u -> dz=-du/u^2 -> estremi:(0,1/20)
# calcolo integrale da 0 a 1/20 di m(y)*g(y) con g() densità di una uniforme (0,1/20) (i.e. g(y)=20 per ogni y in (0,1/20))
# m(y) = 
m <- function(u) {exp(-1/(2*u^2))/(20*(u^2)*sqrt(2*pi))}
g <- function(u) {ifelse(u<=1/20, 20, 0)}

prob.esatta <- integrate(dnorm, lower=20, upper=Inf)$value
y <- runif(R, min=0, max=1/20)
my <- m(y)
stima.mc <- mean(my)
var.mc <- sum((m(y)-stima.mc)^2)/(R*(R-1))
c(prob.esatta, stima.mc)
## [1] 2.753624e-89 2.754502e-89

Metodi per ridurre la varianza dello stimatore MC

Analizziamo 2 metodi per ridurre la varianza dello stimatore MC, pari (per def.) alla media (sugli R valori) della vera varianza della trasformata m(Y).
La stima della varianza dello stimatore MC, come abbiamo già visto, è pari alla media della varianza campionaria. Lo stimatore MC è non distorto e consistente, ma NON è efficiente: esistono altri stimatori non distorti con varianza minore.

METODO 1 : variabili DI CONTROLLO L’idea è quella di definire una o più variabili di controllo n(Y) (di solito i primi k-momenti della v.c. Y) che siano correlate con la trasformata m(Y) e di cui conosciamo il valore atteso.
Calcoliamo c* = -cov(m(Y),n(Y))/var(n(Y)) se non conosciamo a priori la covarianza; quindi applichiamo la formula per lo stimatore MC con variabile di controllo n(Y) :
psi_controllo = psi + c*[n(Y) - mu] ; con mu = E[n(Y)].

set.seed(440)
m <- function(y) {exp(-y)/(1+y^2)}
R <- 10^4
# integrale analitico
integrale <- integrate(m, lower=0, upper=1)$value
# integrale simulato con MC standard
y <- runif(R)
my <- m(y)
stima.mc <- mean(my)
# integrale simulato con MC + variabile di controllo
n <- (1 - y)
mu <- 0.5
var_n <- var(n)
cstar <- -cov(my,n)/var_n
stima.MC.controllo <- mean(stima.mc + cstar*(n - mu))
# confronto
c(integrale, stima.mc, stima.MC.controllo)
## [1] 0.5247971 0.5227909 0.5247129
# riduzione della varianza
var.mc <- var(my)
var.MC.controllo <- var(my + cstar*(n-mu))
c(var.mc, var.MC.controllo)
## [1] 0.059787415 0.001110663
riduzione_relativa <- (1 - var.MC.controllo/var.mc)
riduzione_relativa*100
## [1] 98.14231

Si dimostra che il metodo 1 (variabili di controllo) è analogo ad una regressione di m(Y) sulla variabile di controllo n(Y): la media dei valori fittati corrisponde alla stima.MC.controllo, la varianza stimata dell’errore del modello corrisponde ad R-volte la varianza dello stimatore, mentre l’R-quadro rappresenta la riduzione relativa di varianza rispetto al caso standard.

# equivalenza con il modello di Regressione
set.seed(440)
m <- function(y) {exp(-y)/(1+y^2)}
R <- 10^4
y <- runif(R)
mu <- 0.5
my <- m(y)
ny <- 1 - y

mod <- lm(my ~ ny)
cstar <- -1*mod$coef[2]; cstar
##         ny 
## -0.8400574
stima <- sum(mod$coef*c(1,mu)); stima
## [1] 0.5247129
Rvarianza <- summary(mod)$sigma^2; Rvarianza
## [1] 0.001110774
riduzione_varianza <- summary(mod)$r.squared; riduzione_varianza
## [1] 0.9814231

Stesso metodo, ma visto in un esempio con 2 variabili di controllo. Applichiamo il metodo in entrambi i modi: sia definendo lo stimatore MC modificato con le variabili di controllo, calcolando manualmente i due valori c1, c2 e prendendo la media dei valori; sia come problema di regressione di m(Y) sulla coppia di regressori (variabili di controllo).