Última atualização: 08 de julho, 2025 às 21:39.

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. \[\mathbb{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) \\ \mathbb{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 \(\mathbb{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, \[ \mathbb{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{\mathbb{E}}_p(Y)\) converge para \(\mathbb{E}_p(Y)\) com \(N \rightarrow \infty\) (lei dos grandes números) e possui variabilidade \({\rm Var}[\hat{\mathbb{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*} \mathbb{E}_f[Y] &= \sum_j y_j f(y_j) \approx \frac{1}{N} \sum_i y_i \\ \mathbb{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.499

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.3106

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*} \mathbb{E}_p[\phi(Y)]) &= \sum_j \phi(y_j) p(y_j) \approx \frac{1}{N} \sum_i \phi(y_i)\\ {\rm Var}[\hat{\mathbb{E}}_p[\phi(Y)] &= {\rm Var}[\phi(Y)]/N \approx {\rm Var}[\phi(y_i)]/N. \end{align*}\]

Exemplo 2

Considere a função \(g(y)\) do Exemplo 1 e obtenha a aproximação de Monte Carlo para \(\mathbb{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*} \mathbb{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.8802

Aproximação de esperança pela 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 \(\mathbb{E}_p(Y)\).

A solução de “importance sampling” é: \[\begin{align*} \mathbb{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) = \mathbb{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 é:

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

A média calculada em 4. é uma aproximação para \(\mathbb{E}_p(Y)\).

Em resumo: \[ \mathbb{E}_p(Y) \approx \hat{\mathbb{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 propostas) ou importance distribution (distribuição de importância), instrumental distribution (distribuição instrumental) ou mesmo biased distribution (distribuição com viés).

Exemplo 3

Considere que deseja-se obter a estimativa de \(\mathbb{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.298971

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

A condição requer que

  1. \(p(y)/q(y)\) seja computável,
  2. 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{\mathbb{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. \(\mathbb{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*} \mathbb{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 \[\mathbb{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.321224

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)

sum(yj * py)                ## teórica
## [1] 1.933333
ysim.q <- sample(yj, N, replace = TRUE, prob = qy)
(pqy <- py/qy)              ## pesos
## [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.92598
par(mfrow=c(1,2))
plot(yj, py, type = "h")
plot(yj, qy, type = "h")

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

(TEO <- sum(yj * py))       ## teórica
## [1] 1.933333
ysim.q <- sample(yj, N, replace = TRUE, prob = qy)
(pqy <- py/qy)              ## pesos
## [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.94537
par(mfrow=c(1,2))
plot(yj, py, type = "h")
plot(yj, qy, type = "h")

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

Comparando a eficiê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] 3 3 5 5 3 1 1 3 1 1 ...
##  $ IS1: num [1:10000] 3 3 1 3 1.6 3 3 2.4 1.6 1.6 ...
##  $ IS2: num [1:10000] 0.35 0.35 1.4 0.35 0.35 0.7 0.7 0.7 0.7 0.35 ...
##  $ IS3: num [1:10000] 2.1 1.75 1.75 1.75 1.75 2.1 2.1 2.1 1.75 2.1 ...

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.926400 1.924280 1.956850 1.937157
sapply(sim, function(x) sqrt(var(x)/N))
##          MC         IS1         IS2         IS3 
## 0.011895901 0.009355565 0.024372209 0.001630638

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) \\ \mathbb{E}_p[Y] &= \int y p(f) dy \end{align*}\] O subscrito na notação \(\mathbb{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, \[ \mathbb{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{\mathbb{E}}_p(Y)\) converge para \(\mathbb{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*} \mathbb{E}_p(Y) &= \int y p(y) dy = \int \left(y \frac{p(y)}{q(y)}\right) q(y) = \mathbb{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:

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

A média calculada em 4. é uma aproximação para \(\mathbb{E}_p(Y)\).

Em resumo: \[ \mathbb{E}_p(Y) \approx \hat{\mathbb{E}}_q\left(Y \frac{p(Y)}{q(Y)}\right). \] De forma mais geral, \[ \mathbb{E}_p(\phi(Y)) = \int \phi(y) p(y) dy \approx \hat{\mathbb{E}}_q\left(\phi(Y) \frac{p(Y)}{q(Y)}\right) = \frac{1}{N} \sum_i \phi(y_i) w_i .\] e o 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 \(\mathbb{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] 3.005841
mean(samQ1*w1)   ## amostragem por importancia com proposal uniforme
## [1] 2.970598
mean(samQ2*w2)   ## amostragem por importancia com proposal normal
## [1] 3.006236

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.00817665
sqrt(var(samQ1*w1)/N)   ## amostragem por importancia com proposal uniforme
## [1] 0.03120326
sqrt(var(samQ2*w2)/N)   ## amostragem por importancia com proposal normal
## [1] 0.01123742

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.

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

  1. Obter \(\mathbb{E}_p[Y]\) com \(p(Y) \sim \text{Beta}(3,5) \text{ e } q(Y) \sim {\rm Beta}(6,2)\).

  2. Obter \(\mathbb{E}_p[Y]\) com \(p(Y) \sim {\rm HalfNormal}(1) \text{ e } q(Y) \sim {\rm Exp}(2)\).

  3. Obter \(P_p[Y > 4.5]\) com \(p(Y) \sim {\rm N}(0,1) \text{ e } q(Y) \sim {\rm ShiftedExp}(1, 4.5)\).

  4. Deseja-se obter \(\mathbb{E}_p[\exp\{-Y + {\rm 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 de Monte Carlo.
  • O que se obtém usando agora \(q(Y) \sim {\rm Exp}(2,5)\)?

OBS: Expressões das densidades de probabilidades

  • \(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\)

Licença Creative Commons 4.0

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