Ce TD illustre une technique assez bluffante pour réduire la variance et donc la qualité de ses estimations. Nous allons l'illustrer en sur un problème simple, l'estimation de \( \pi \) par la méthode de Monte-Carlo. On tire un point uniformément dans le carré unitaire et s'il est à l'intérieur du quart de cercle, on compte 1 et 0 sinon. La variable aléatoire suit une loi de Bernouilli de paramètre \( \pi/4 \) et est donc d'espérance \( \frac{\pi}{4} \) et de variance \( \frac{\pi}{4}(1-\frac{\pi}{4}) \).
Si vous réalisez la moyenne de \( n \) tirages, quelle est l'intervalle de confiance associé à cet estimateur?
library(plyr)
library(ggplot2)
pi_estim = function(n = 100) {
X = runif(n)
Y = runif(n)
T = X^2 + Y^2
U = ifelse(T <= 1, 1, 0)
data.frame(X = X, Y = Y, U = U, T = T)
}
df = pi_estim(100)
head(df)
## X Y U T
## 1 0.03588 0.10188 1 0.01167
## 2 0.71957 0.95075 0 1.42171
## 3 0.65677 0.01082 1 0.43146
## 4 0.66844 0.36293 1 0.57852
## 5 0.86396 0.64145 0 1.15789
## 6 0.47337 0.04959 1 0.22654
ggplot(data = df, aes(x = X, y = Y, color = factor(U))) + geom_point() + coord_fixed() +
theme_bw()
ggplot(data = df, aes(x = X, y = U)) + geom_point()
Erm, normalement, c'est comme ça qu'il faut faire quand on veut voir la correlation entre deux variables mais là ça n'est pas très lisible car \( Y \) est dans \( \{0,1\} \). Essayons avec un histogramme pour mieux voir comment ça se réparti:
ggplot(data = df, aes(x = X)) + geom_histogram() + facet_wrap(~U, ncol = 1)
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
Là, du coup, on voit bien que quand \( X \) est proche de \( 0 \), il y a très peu de chances pour que \( U \) vaille 0 (i.e., que le point ne soit pas dans le disque), ce qui est assez logique.
ggplot(data = df, aes(x = X^2, y = U)) + geom_point()
ggplot(data = df, aes(x = X^2)) + geom_histogram() + facet_wrap(~U, ncol = 1)
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
Mouaip, bof.
ggplot(data = df, aes(x = X^2 + Y^2, y = U)) + geom_point()
ggplot(data = df, aes(x = X^2 + Y^2)) + geom_histogram() + facet_wrap(~U, ncol = 1)
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
Qu'en déduisez-vous ?
On s'intéresse à la variable \( V=U+\alpha.(X-\frac{1}{2}) \) (pour un \( \alpha \) fixé).
On peut montrer sans grande difficulté que \( var(V)=Var(X).\alpha^2 + 2cov(X,U).\alpha + var(U) \). En choisissant astucieusement \( \alpha \), il est donc possible de diminuer \( var(V) \) et donc d'avoir une meilleur estimation à partir d'exactement les mêmes informations!!!
df = pi_estim(100)
var(df$X)
## [1] 0.07402
cov(df$X, df$U)
## [1] -0.06566
var(df$U)
## [1] 0.1991
alpha_est = -cov(df$X, df$U)/var(df$X)
alpha_est
## [1] 0.887
N = 100
df = data.frame()
for (i in (1:N)) {
d = pi_estim(100)
d$trial = i
for (a in c(0, 0.7)) {
d$alpha = a
d$V = d$U + a * (d$X - 1/2)
df = rbind(df, d)
}
}
df_ci = ddply(df, c("trial", "alpha"), summarise, est = mean(V), var_est = var(V),
conf = 2 * sd(V)/sqrt(length(V)))
ggplot(data = df_ci, aes(x = trial, y = var_est, color = factor(alpha))) + geom_point() +
ylim(0, max(df_ci$var_est))
ggplot(data = df_ci, aes(x = trial, y = est, ymin = est - conf, ymax = est +
conf, color = factor(alpha))) + geom_point() + geom_errorbar() + ylim(0.5,
1) + geom_hline(yintercept = pi/4)
# + facet_wrap(~factor(alpha))
df = pi_estim(100)
var(df$T)
## [1] 0.1913
cov(df$T, df$U)
## [1] -0.1214
var(df$U)
## [1] 0.1425
alpha_est = -cov(df$T, df$U)/var(df$T)
alpha_est
## [1] 0.6349
N = 100
df = data.frame()
for (i in (1:N)) {
d = pi_estim(100)
d$trial = i
# crude
d$method = "crude"
d$V = d$U
df = rbind(df, d)
# smarter
d$method = "smarter"
d$V = d$U + 0.7 * (d$X - 1/2)
df = rbind(df, d)
# super smart
d$method = "super smart"
d$V = d$U + 0.7 * (d$T - 2/3)
df = rbind(df, d)
}
df_ci = ddply(df, c("trial", "method"), summarise, est = mean(V), var_est = var(V),
conf = 2 * sd(V)/sqrt(length(V)))
ggplot(data = df_ci, aes(x = trial, y = est, ymin = est - conf, ymax = est +
conf, color = method)) + geom_point() + geom_errorbar() + geom_hline(yintercept = pi/4) +
facet_wrap(~method)
Toujours pas convaincus ? Regardons sur l'ensemble des échantillons alors pour que ça soit plus lisible.
df_ci = ddply(df, c("method"), summarise, est = mean(V), var_est = var(V), conf = 2 *
sd(V)/sqrt(length(V)))
ggplot(data = df_ci, aes(x = method, y = est, ymin = est - conf, ymax = est +
conf, color = method)) + geom_point() + geom_errorbar() + geom_hline(yintercept = pi/4)
Mettons que l'on souhaite estimer \( \int_0^1 \frac{1}{1+x}dx \) (comme sur ). On peut regarder l'espérance de \( \frac{1}{1+X} \) avec \( X \) uniforme dans \( [0,1] \).
x = runif(100)
y = 1/(1 + x)
mean(y)
## [1] 0.6833
var(y)
## [1] 0.02081
Évidemment la variance est relativement élevée et l'intervalle de confiance du coup relativement large. On notera que 1-x est également un échantillon tiré uniformément dans \( [0,1] \) et que si jamais on a “pas eu de chance” et qu'on a tiré trop de points vers le 0, 1-x aurait la propriété inverse, ce qui rééquilibrerait un peu les choses pour la moyenne. Regardons donc \( Y=\frac{1}{2}\left(\frac{1}{1+X}+\frac{1}{1+1-X}\right) \). \( Y \) a évidemment la même espérance que \( X \), mais une variance bien bien plus faible!
y = 1/2 * (1/(1 + x) + 1/(1 + 1 - x))
mean(y)
## [1] 0.6945
var(y)
## [1] 0.0007068
Appliquons ça à notre calcul de \( \pi \):
pi_estim_antithetic = function(n = 100) {
d = data.frame()
X = runif(n)
Y = runif(n)
T = X^2 + Y^2
U = ifelse(T <= 1, 1, 0)
X = 1 - X
T = X^2 + Y^2
U = U + ifelse(T <= 1, 1, 0)
Y = 1 - Y
T = X^2 + Y^2
U = U + ifelse(T <= 1, 1, 0)
X = 1 - X
T = X^2 + Y^2
U = U + ifelse(T <= 1, 1, 0)
U = U/4
d = rbind(d, data.frame(X = X, Y = Y, U = U, T = T))
d
}
Résultat:
df = pi_estim(100)
mean(df$U)
## [1] 0.81
var(df$U)
## [1] 0.1555
df2 = pi_estim_antithetic(100)
mean(df2$U)
## [1] 0.795
var(df2$U)
## [1] 0.02826
En simulation, évidemment, ça vaut vraiment le coup mais on peut éventuellement utiliser le même principe pour des mesures expérimentales, Pas simple mais faisable…