Última atualização: 06 de junho, 2025 às 18:19.

Amostragem por importância (importance sampling) é um método de Monte Carlo para avaliar propriedades de uma distribuição a partir de amostras obtidas de outra distribuição. A motivação pode ser

obter resultados quando e amostragem da distribuição original é difícil ou impossível,
reduzir o erro de Monte Carlo.

Variável aleatória discreta

Aproximação de esperança por Monte Carlo

Objetivo: Obter o valor esperado de uma v.a. $Y$ discreta. \[E[Y] = \sum_j y_j P[Y = y_j] \]

Vamos adotar seguinte notação para distribuição da variável e seu valor esperado: \[\begin{align*} Y &\sim p(y) \\ E_p[Y] &= \sum_j y_j P[Y = y_j] = \sum_j y_j p(y_j). \end{align*}\] O subscrito na notação $E_p[Y]$ destaca que o valor esperado é sob a distribuição $p(y)$.

Este valor esperado pode ser aproximado por simulação de Monte Carlo, \[ E_p[Y] = \sum_j y_j p(y_j) \approx \frac{1}{N} \sum_i y_i ,\] ou seja, a média de amostras tomadas da distribuição $p(y)$. O subscrito $j$ em $y_j$ se refere as possíveis valores que a variável pode assumir, enquanto que o subscrito $j$ em $y_j$ indica um valor dentre o $N$ simulados da distribuição.

A estimativa $\hat{E}_p(Y)$ converge para $E_p(Y)$ com $N \rightarrow \infty$ (lei dos grandes números) e possui variabilidade ${\rm Var}[\hat{E}_p(Y)] = {\rm Var}[Y]/N$.

Exemplo 1

Sejam duas funções de probabilidades discretas com gráficos dados a seguir. O objetivo é obter uma aproximação de Monte Carlo da esperança da variável em relação a cada uma das duas distribuições dadas. \[\begin{align*} f(y_j) &= \frac{1}{6} \;\; \forall j \in \{1, 2, 3, 4, 5, 6\} \\ g(y_j) &= \frac{i}{21} \;\; \forall j \in \{1, 2, 3, 4, 5, 6\} \end{align*}\]

yj <- 1:6
fy <- rep(1/6, 6)
gy <- yj/21
par(mfrow=c(1,2), mar = c(3,3,0.5, 0.5), mgp = c(1.8, 0.8, 0))
plot(yj, fy, type = "h")
plot(yj, gy, type = "h")

As esperanças exatas e aproximadas numericamente são \[\begin{align*} E_f[Y] &= \sum_j y_j f(y_j) \approx \frac{1}{N} \sum_i y_i \\ E_g[Y] &= \sum_j y_j g(y_j) \approx \frac{1}{N} \sum_i y_i \end{align*}\]

N <- 1e4              ## número de simulações 

sum(yj * fy)          ## teórica
## [1] 3.5
ysim.f <- sample(yj, N, replace = TRUE, prob = fy)
mean(ysim.f)          ## aproximação por simulação
## [1] 3.5137

sum(yj * gy)          ## teórica
## [1] 4.333333
ysim.g <- sample(yj, N, replace = TRUE, prob = gy)
mean(ysim.g)          ## aproximação por simulação de g
## [1] 4.3501

O método é válido não apenas para aproximar o valor esperado de $Y$ mas também para alguma função $\phi(Y)$ e temos portanto \[\begin{align*} E_p[\phi(Y)]) &= \sum_j \phi(y_j) p(y_j) \approx \frac{1}{N} \sum_i \phi(y_i)\\ {\rm Var}[\hat{E}_p[\phi(Y)] =& {\rm Var}[\phi(Y)]/N. \end{align*}\]

Exemplo 2

Considere a função $g(y)$ do Exemplo 1 e obtenha a aproximação de Monte Carlo para $E_g[Y^2]$, ou seja, o valor esperado de $\phi(Y) = Y^2$ com $N = 10000$ simulações.

A esperança exata e aproximada numericamente são \[\begin{align*} E_g[Y^2] &= \sum_j y^2_j g(y_j) \approx \frac{1}{N} \sum_i y_i^2 \end{align*}\]

N <- 1e4              ## número de simulações 
yj <- 1:6
gy <- yj/21

sum(yj^2 * gy)         ## teórica
## [1] 21

ysim.g <- sample(yj, N, replace = TRUE, prob = gy)
mean(ysim.g^2)         ## aproximação por simulação de g
## [1] 20.8236

Aproximação de esperança por amostragem por importância

A aproximação de Monte Carlo requer que seja possível tomar amostras de $p(y)$. Além disto, deseja-se que o erro de Monte Carlo seja a menor possível. O método da amostragem por importância consiste em obter amostras de uma outra distribuição, $q(y)$, que são devidamente ponderadas para, ainda assim, obter a estimativa de $E_p(Y)$.

A solução de “importance sampling” é: \[\begin{align*} E_p(Y) &= \sum_j y_j p(y_j) = \sum_j \left(y_j \frac{p(y_j)}{q(y_j)}\right) q(y_j) = E_q\left( y_j \frac{p(y_j)}{q(y_j)} \right) \\ &\approx \frac{1}{N} \sum_i y_i \frac{p(y_i)}{q(y_i)} = \frac{1}{N} \sum_i y_i w_i , \end{align*}\] em que os $y_i$ são amostras de $q(y)$. Portanto, suponha que não se sabe como amostrar de $p(y)$ ou deseja-se reduzir erro Monte Carlo e sabemos como amostrar de uma outra distribuição $q(y)$. Podemos calcular para cada valor $y_i$ amostrado $q(y)$ a quantidade $w_i = p(y_i)/q(y_i)$, chamada de ponderação (peso) de importância. Em resumo, o algoritmo é:

simule $y_i$ de $q(y_i)$,
calcule $y_i \frac{p(y_i)}{q(x_i)} = y_i w_i$,
repita 1-2 muitas vezes,
Calcule a média dos valores calculados em 2.

A média calculada em 4. é uma aproximação para $E_p(Y)$.

Em resumo: \[ E_p(Y) \approx \hat{E}_q\left(Y \frac{p(Y)}{q(Y)}\right). \]

A função $q(y)$ da qual tomamos amostras é chamada de proposal distribution (distribuição de proposta) ou importance distribution (distribuição de importância) ou mesmo biased distribution (distribuição com viés).

Exemplo 3

Considere que deseja-se obter a estimativa de $E_p(Y)$, porém a partir de amostras geradas de $q(y)$.

py <- yj/21
qy <- rep(1/6, 6)
ysim.q <- sample(yj, N, replace = TRUE, prob = qy)
(pqy <- py/qy)
## [1] 0.2857143 0.5714286 0.8571429 1.1428571 1.4285714 1.7142857
mean(ysim.q * pqy[ysim.q])  # aproximação por importance sampling
## [1] 4.284914

O resultado pode ser comparado com valor teórico 1.9333333.

A condição requer que

$p(y)/q(y)$ seja computável,
o domínio de $p(y)$ esteja contido no domínio de $q(y)$.

Note que $p(y)/q(y)$ são pesos aplicados aos valores simulados e portanto $\hat{E}_p(Y)$ é uma média ponderada de valores simulados de $q(y)$.

Generalizações:

aplica-se não apenas para valor esperado mas também para funções da v.a. $E_p(\phi(y))$,
aplica-se de forma análoga para distribuições contínuas.

Observação

A função objetivo $p(y)$ pode ser conhecida “a uma constante de proporcionalidade” ou seja, $p(y) = k(y_i)/C$ e ainda assim a amostragem por importância se aplica. Considere que \[\begin{align*} p(y_j) &\propto j \;\; \forall j \in {1, 2, 3, 4, 5, 6}\\ q(y_j) &= \frac{1}{6} \;\; \forall j \in {1, 2, 3, 4, 5, 6} \end{align*}\]

Tem-se que \[\begin{align*} E_p[Y] &= \sum_j y_j p(y_j) = \sum_j y_j \frac{k(y_j)}{C} = \sum_j y_j \frac{k(y_j)}{C q(y_j)} q(y_j) \\ &\approx \frac{1}{C}\frac{1}{N} \sum_i y_i \frac{k(y_i)}{q(y_i)} = \frac{1}{C}\frac{1}{N} \sum_i y_i w_i , \end{align*}\] o que remete a como calcular a contante $C$, que é simplesmente e média dos pesos $\overline{w} = \sum_i w_i/N$ pois \[\begin{align*} \sum_j p(y_j) = \sum_j \frac{k(y_j)}{C} &= 1 \\ C &= \sum_j k(y_j) = \sum_j \frac{k(y_j)}{q(y_j)} q(y_j)\\ & \approx \frac{1}{N} \sum_i \frac{k(y_i)}{q(y_i)} = \sum_i w_i = \overline{w} \end{align*}\] Portanto retornando à expressão anterior \[E_p[Y] \approx \frac{1}{C}\frac{1}{N} \sum_i y_i w_i = \frac{1}{\overline{w}}\frac{1}{N} \sum_i y_i w_i = \frac{\overline{y w}}{\overline{w}}\]

O código a seguir ilustra os passos de amostragem por importância nesta situação.

N <- 1e4
yj <- 1:6

ky <- 1:6
qy <- rep(1/6, 6)

ysim.q <- sample(yj, N, replace = TRUE, prob = qy)
(kqy <- ky/qy)
## [1]  6 12 18 24 30 36
wi <- kqy[ysim.q]  
mean(ysim.q*wi)/mean(wi)
## [1] 4.335906

Exemplo 4

Mudando agora a função objetivo $p(y)$

$y_j$	1	2	3	4	5	6
$p(y_j)$	1/2	1/4	2/15	1/15	1/30	1/60
$q(y_j)$	1/6	1/6	1/6	1/6	1/6	1/6
$w_j = p(y_j)/q(y_j)$	3	1.5	0.8	0.4	0.2	0.1

N <- 1e4     
yj <- 1:6

py <- c(1/2, 1/4, 2/15, 1/15, 1/30, 1/60)
qy <- rep(1/6, 6)

par(mfrow=c(1,2))
plot(yj, py, type = "h")
plot(yj, qy, type = "h")

sum(yj * py)                ## teórica

## [1] 1.933333

ysim.q <- sample(yj, N, replace = TRUE, prob = qy)
(pqy <- py/qy)

## [1] 3.0 1.5 0.8 0.4 0.2 0.1

mean(ysim.q * pqy[ysim.q])  ## aproximação por importance sampling

## [1] 1.9318

Exemplo 5

Mudando agora a função de importância $q(y)$

$y_j$	1	2	3	4	5	6
$p(y_j)$	1/2	1/4	2/15	1/15	1/30	1/60
$q(y_j)$	1/21	2/21	3/21	4/21	5/21	6/21
$w_j = p(y_j)/q(y_j)$	3	1.5	0.8	0.4	0.2	0.1

N <- 1e4     
yj <- 1:6

py <- c(1/2, 1/4, 2/15, 1/15, 1/30, 1/60)
qy <- yj/21

par(mfrow=c(1,2))
plot(yj, py, type = "h")
plot(yj, qy, type = "h")

(TEO <- sum(yj * py))                ## teórica

## [1] 1.933333

ysim.q <- sample(yj, N, replace = TRUE, prob = qy)
(pqy <- py/qy)

## [1] 10.50000000  2.62500000  0.93333333  0.35000000  0.14000000  0.05833333

mean(ysim.q * pqy[ysim.q])  ## aproximação por importance sampling

## [1] 1.924335

Neste exemplo a aproximação pode parecer ineficiente mas ainda funciona!

Comparando a eficência de opções para $q(y)$

Vamos fazer um gráfico da função objetivo e três opções para função de importancia.

py <- c(1/2, 1/4, 2/15, 1/15, 1/30, 1/60)
q1y <- rep(1/6, 6)
q2y <- (1:6)/21
q3y <- (6:1)/21
par(mfrow=c(1,3))
plot((1:6)+0.1, py, type = "h", xlab = "y", xlim = c(0.5, 6.5))
lines((1:6)-0.1, q1y, type = "h", col = "blue")
plot((1:6)+0.1, py, type = "h", xlab = "y", xlim = c(0.5, 6.5))
lines((1:6)-0.1, q2y, type = "h", col = "blue")
plot((1:6)+0.1, py, type = "h", xlab = "y", xlim = c(0.5, 6.5))
lines((1:6)-0.1, q3y, type = "h", col = "blue")

Vamos agora fazer a aproximação Monte Carlo e as três de amostragem por importância. No código em R definir uma função para facilitar a amostragem poor importância. Criamos uma lista para armazenar as simulações de cada método.

funIS <- function(q, p, yvals = 1:6, N = 1e4, retorna = c("media", "sims")){
    retorna <- match.arg(retorna)
    ysim <- sample(yvals, N, replace = TRUE, prob = q)
    pq <- p/q
    yw <- ysim * pq[ysim]
    if(retorna == "sims") return(yw)
    else return(mean(yw)) 
}
sim <- list()
sim$MC <- sample(yj, N, prob = py, replace = TRUE)
sim$IS1 <- funIS(q = q1y, p = py, retorna = "sims")
sim$IS2 <- funIS(q = q2y, p = py, retorna = "sims")
sim$IS3 <- funIS(q = q3y, p = py, retorna = "sims")
str(sim)
## List of 4
##  $ MC : int [1:10000] 1 5 4 4 1 1 3 1 1 3 ...
##  $ IS1: num [1:10000] 3 1.6 3 2.4 0.6 3 1.6 0.6 1.6 0.6 ...
##  $ IS2: num [1:10000] 1.4 0.7 0.7 0.35 0.35 1.4 0.7 0.35 0.35 0.35 ...
##  $ IS3: num [1:10000] 1.75 1.75 1.75 2.1 1.75 ...

Ao final extraímos as estimativas e os respectivos erros padrão, o que permite comparar a eficiências de cada uma das aproximações.

sapply(sim, mean)
##       MC      IS1      IS2      IS3 
## 1.928500 1.932800 1.903230 1.932992
sapply(sim, function(x) sqrt(var(x)/N))
##          MC         IS1         IS2         IS3 
## 0.012090219 0.009374882 0.023860601 0.001633460

Lembrando que o valor teórico da esperança é 1.9333.

Variável aleatória contínua

De forma análoga ao caso discreto, o valor esperado de uma v.a. $Y$ contínua é \[\begin{align*} Y &\sim p(y) \\ E_p[Y] &= \int y p(f) dy \end{align*}\] O subscrito na notação $E_p[Y]$ destaca que o valor esperado é sob a distribuição $p(y)$.

Este valor esperado pode ser aproximado por simulação de Monte Carlo, \[ E_p(Y) = \sum_j y_i p(y_j) \approx \frac{1}{N} \sum_i y_i ,\] ou seja, a média de amostras tomadas da distribuição $p(y)$.

A estimativa $\hat{E}_p(Y)$ converge para $E_p(Y)$ com $N \rightarrow \infty$ (lei dos grandes números) e possui variabilidade ${\rm Var}[\hat{E}_p(Y)] = {\rm Var}[Y]/N$.

A solução de “importance sampling” fica: \[\begin{align*} E_p(Y) &= \int y p(y) dy = \int \left(y \frac{p(y)}{q(y)}\right) q(y) = E_q\left(y \frac{p(y)}{q(y)} \right) \\ &\approx \frac{1}{N} \sum_i y_i \frac{p(y_i)}{q(y_i)} = \frac{1}{N} \sum_i y_i w_i , \end{align*}\] em que os $y_i$ são amostras de $q(y)$. Portanto, suponha que não se sabe como amostrar de $p(y)$ ou deseja-se reduzir a variância de Monte Carlo, mas sabemos como amostrar de $q(y)$. Podemos calcular para cada valor $y_i$ amostrado $q(y)$ a quantidade $w_i = p(y_i)/q(y_i)$, chamada de ponderação (peso) de importância. O algoritmo é o mesmo do caso discreto:

simule $y_i$ de $q(y_i)$,
calcule $w_i \frac{p(y_i)}{q(x_i)}$,
repita 1-2 muitas vezes,
Calcule a média ponderada $y_i w_i$ dos valores calculados.

A média calculada em 4. é uma aproximação para $E_p(Y)$.

Em resumo: \[ E_p(Y) \approx \hat{E}_q\left(Y \frac{p(Y)}{q(Y)}\right). \] De forma mais geral, \[ E_p(\phi(Y)) = \int \phi(y) p(f) dy \approx \hat{E}_q\left(\phi(Y) \frac{p(Y)}{q(Y)}\right) = \frac{1}{N} \sum_i \phi(y_i) w_i .\] e a erro padrão da aproximação (estimativa) é \[ \sqrt{\text{Var}[\hat{E}_p(\phi(Y))]} = \sqrt{\frac{\text{Var}[\phi(Y) w]}{N}}. \]

Exemplo 6

Considere que deseja-se obter aproximar por simulação $E_p(Y)$ através de importance sampling considerando duas distribuições de importãncia.

$p(y)$ é dada por $Y \sim Tri(a=1, b = 5, c = 3)$
$q_1(y)$ é dada por $Y \sim U(0, 6)$
$q_2(y)$ é dada por $Y \sim N(3)$

N <- 1e4
samP <- rtriang(N, a = 1, b = 5, c = 3)
samQ1 <- runif(N, 0, 6)
samQ2 <- rnorm(N, 3, 1)

w1 <- dtriang(samQ1, a = 1, b = 5, c = 3)/dunif(samQ1, 0, 6)
w2 <- dtriang(samQ2, a = 1, b = 5, c = 3)/dnorm(samQ2, 3, 1)

mean(samP)       ## aproximação de Monte Carlo
## [1] 2.993339
mean(samQ1*w1)   ## amostragem por importancia com proposal uniforme
## [1] 3.021907
mean(samQ2*w2)   ## amostragem por importancia com proposal normal
## [1] 2.992835

Pode-se comparar os erros padrão das aproximações. Neste exemplo, as amostragem por importância não melhoraram o erro padrão da aproximação. Entre as duas opções para função de importância nota-se que a segunda fornece estimativa com menor erro padrão.

sqrt(var(samP)/N)       ## aproximação de Monte Carlo

## [1] 0.008185623

sqrt(var(samQ1*w1)/N)   ## amostragem por importancia com proposal uniforme

## [1] 0.03131184

sqrt(var(samQ2*w2)/N)   ## amostragem por importancia com proposal normal

## [1] 0.01132287

Agora experimente com variações deste exemplo e proponha e implemente outras opções de $p(y)$ e $q(y)$!

Exercícios

Nos exercícios 1-3 a seguir siga os seguintes itens.

Obtenha o resultado teórico de cada quantidade de interessse.
Escreva um algoritmo que pode ser usado para aproximar $E_p[\phi(Y)]$ por simulação de Monte Carlo.
Escreva um algoritmo que pode ser usado para aproximar $E_p[\phi(Y)]$ por amostragem por importância.
Inclua no código gráficos das funções envolvidas.
Discuta a intuição, resultados e interpretação em cada caso.

Obter $E_p[Y]$ com $p(Y) \sim \text{Beta}(3,5) \text{ e } q(Y) \sim {\rm Beta}(6,2)$.
Obter $E_p[Y]$ com $p(Y) \sim {\rm HalfNormal}(1) \text{ e } q(Y) \sim {\rm Exp}(2)$.
Obter $P_p[Y > 4.5]$ com $p(Y) \sim {\rm N}(0,1) \text{ e } q(Y) \sim {\rm ShiftedExp}(1, 4.5)$.
Deseja-se obter $E_p[\exp\{-Y + cos(Y)\}]$ com $p(Y) \sim {\rm Exp}(1) \text{ e } q(Y) \sim {\rm Exp}(2)$.
- Obter o resultado teórico por algum método de aproximação numérica (sugestão: usar função integrate() do R).
- Obter a aproximação por Monte Carlo.
- Obter a aproximação por importance sampling.
- Obtenha o erro Monte Carlo de cada método.
- Faça um gráfico com as funções $\phi(y)p(y)$, $p(y) e $q(y)$ e associe comportamentos das funções com os erros Monte Carlo.

O que se obtém usando agora $q(Y) \sim {\rm Exp}(2,5)$?

OBS: Densidades:

$Y \sim \text{Beta}(\alpha, \beta) \longrightarrow f(y) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} y^{\alpha - 1} (1 - y)^{\beta - 1}, 0 < y < 1$
$Y \sim \text{Normal}(\mu,\sigma) \longrightarrow f(y) = \frac{1}{\sigma \sqrt{2\pi}} \exp\left( -\frac{(y-\mu)^2}{2\sigma^2} \right), -\infty < y < \infty$
$Y \sim \text{Half-Normal}(\sigma) \longrightarrow f(y) = \frac{\sqrt{2}}{\sigma \sqrt{\pi}} \exp\left( -\frac{y^2}{2\sigma^2} \right), y \geq 0$
$Y \sim {\text{Exp}}(\lambda) \longrightarrow f(y) = \lambda e^{-\lambda y}, y \geq 0$
$Y \sim {\text{Shifted Exp}}(\lambda, \mu) \longrightarrow f(y) = \lambda e^{-\lambda (y - \mu)}, y \geq \mu$

Este conteúdo está disponível por meio da Licença Creative Commons 4.0

Amostragem por importância (importance sampling)

Paulo Justiniano