Amostragem preferencial

A amostragem preferencial na amostragem de importância em inglês é um método de redução de variância que pode ser usado no método de Monte Carlo . A ideia subjacente à amostragem preferencial, EP a seguir, é que certos valores tomados por uma variável aleatória em uma simulação têm mais efeito do que outros no estimador desejado. Se esses valores grandes ocorrerem com mais frequência, a variância de nosso estimador pode ser reduzida.

Portanto, o método PE é escolher uma distribuição que “incentive” os valores importantes. Usar uma distribuição enviesada levará a um estimador enviesado se o aplicarmos diretamente às simulações. No entanto, as diferentes simulações são ponderadas a fim de corrigir esse viés; o estimador EP é então imparcial. O peso dado a cada simulação é a razão de verossimilhança, que é a densidade Radon-Nikodym da distribuição verdadeira versus a distribuição enviesada.

O ponto fundamental na implementação de uma simulação utilizando o PE é a escolha da distribuição enviesada. Escolher ou criar um bom elenco tendencioso é a arte dos EPs. A vantagem pode ser uma grande economia de tempo computacional, enquanto a desvantagem de uma má distribuição pode ser cálculos mais longos do que uma simples simulação de Monte-Carlo.

Teoria

Em Monte-carlo

Queremos estimar uma quantidade G , que é expressa na forma de uma integral:

{\ displaystyle G = \ int _ {a} ^ {b} g (x) \, {\ mbox {d}} x}

Consideramos aqui uma integração na dimensão 1, mas podemos generalizar para qualquer dimensão.

O princípio básico dos métodos de Monte-Carlo é ver a integral anterior como

{\ displaystyle G = (ba) \ int _ {a} ^ {b} g (x) f_ {X} (x) \, {\ mbox {d}} x = (ba) \, \ mathbb {E} (g (X))}

onde X é uma variável aleatória uniformemente distribuída em [ a; b ] e sua densidade. ${\ displaystyle f_ {X} (\ cdot) = {\ frac {1} {ba}}}$

Se tivermos uma amostra , distribuída de forma idêntica e independente (iid) de acordo com , podemos estimar G por: ${\ displaystyle (x_ {1}, x_ {2}, \ cdots, x_ {N})}$ ${\ displaystyle {\ mathcal {U}} ([a; b])}$

{\ displaystyle {\ hat {g}} _ {N} = {\ frac {(ba)} {N}} \ sum _ {i = 1} ^ {N} g (x_ {i})}

É um estimador imparcial (isto é ) e consistente (de acordo com a lei dos grandes números) de G. Sua variação é: ${\ displaystyle \ mathbb {E} ({\ hat {g}} _ {N}) = G}$

{\ displaystyle \ sigma _ {{\ hat {g}} _ {N}} ^ {2} = {\ frac {(ba) ^ {2} \ sigma _ {g} ^ {2}} {N}} }

com a variância da variável aleatória ${\ displaystyle \ sigma _ {g} ^ {2}}$ $g (X)$

{\ displaystyle \ sigma _ {g} ^ {2} = {\ frac {1} {(ba)}} \ int _ {a} ^ {b} g ^ {2} (x) \, {\ mbox { d}} x- \ left ({\ frac {1} {ba}} \ int _ {a} ^ {b} g (x) \, {\ mbox {d}} x \ right) ^ {2}}

Princípio da amostragem preferencial

A ideia principal da amostragem preferencial é substituir na simulação a densidade uniforme por uma densidade alternativa (ou densidade enviesada ), notada , que tenta imitar a função g . Com isso, as estampas uniformes, que não favorecem nenhuma região, são substituídas por estampas mais “fiéis”. Assim, a amostragem é feita de acordo com a importância da função g : é inútil desenhar nas regiões onde g assume valores não significativos, para, pelo contrário, concentrar-se nas regiões de alta importância. Esperamos reduzir a variância dessa forma . Em outras palavras, se um dado nível de erro é fixo, a amostragem preferencial permite reduzir teoricamente o número de simulações N em relação a um método clássico de Monte-Carlo. $[a; b]$ ${\ displaystyle f ^ {\ ast} \,}$ ${\ displaystyle \ sigma _ {g} ^ {2}}$

A integral a ser estimada é reescrita como:

{\ displaystyle G = \ int _ {a} ^ {b} {\ frac {g (x)} {f ^ {\ ast} (x)}} f ^ {\ ast} (x) \, {\ mbox {d}} x}

que equivale a:

{\ displaystyle G = \ mathbb {E} ^ {\ ast} [w (X)]}

onde colocamos (chamada de razão de verossimilhança ) e onde X é simulado de acordo com a densidade . É fácil generalizar os resultados anteriores: o estimador de G é ${\ displaystyle w (x) = g (x) / f ^ {\ ast} (x)}$ ${\ displaystyle f ^ {\ ast}}$

{\ displaystyle {\ tilde {g}} _ {N} = {\ frac {1} {N}} \ sum _ {i = 1} ^ {N} w (x_ {i})}

onde é uma amostra iid de acordo com a densidade . A variância do estimador é dada por ${\ displaystyle (x_ {1}, x_ {2}, \ cdots, x_ {N})}$ ${\ displaystyle f ^ {\ ast}}$

{\ displaystyle {\ mbox {Var}} ^ {\ ast} ({\ tilde {g}} _ {N}) = {\ frac {{{\ mbox {Var}} ^ {\ ast} [w (X )]} {NÃO}}}

com finalmente

{\ displaystyle {\ mbox {Var}} ^ {\ ast} [w (X)] = {\ mbox {Var}} ^ {\ ast} \ left [{\ frac {g (X)} {f ^ { \ ast} (X)}} \ right] = \ int _ {a} ^ {b} \ left [{\ frac {g (x)} {f ^ {\ ast} (x)}} \ right] ^ {2} f ^ {\ ast} (x) \, {\ mbox {d}} xG ^ {2}}

Portanto, o problema é se concentrar em obter uma densidade enviesada de forma que a variância do estimador EP seja menor que a do método clássico de Monte-Carlo. A densidade que minimiza a variância (o que a torna zero sob certas condições) é chamada de densidade polarizada ótima . Este último é igual a ${\ displaystyle f ^ {\ ast} \,}$

{\ displaystyle f ^ {\ ast} (x) = {\ frac {| g (x) |} {\ displaystyle \ int _ {a} ^ {b} | g (x) | \, {\ mbox {d }} x}}}

mas esta escolha é ineficaz, porque procuramos precisamente o denominador. No entanto, pode-se esperar reduzir a variância escolhendo uma densidade reproduzindo a função g . ${\ displaystyle f ^ {\ ast}}$

Quase monte carlo

Para estimar a integral , também podemos dispensar todo o formalismo probabilístico anterior. Em vez de usar variáveis aleatórias, usamos sequências de baixa discrição (sequências de Sobol, por exemplo). Na dimensão 1, a abordagem mais simples é ${\ displaystyle G = \ int _ {a} ^ {b} g (x) \, {\ mbox {d}} x}$

{\ displaystyle G = \ int _ {a} ^ {b} g (x) \, {\ mbox {d}} x \ quad \ simeq \ quad {\ frac {ba} {N}} \ sum _ {i = 1} ^ {N} g \ left (a + (ba) {\ frac {i} {N}} \ right)}

Como no Monte Carlo usual, essa aproximação da integral converge tanto mais rápido quanto a função g se aproxima de uma constante. Se g for rigorosamente constante, basta tomar N = 1 para ter a integral exata. Reduzir a variância de g também é crucial aqui; para este fim, a amostragem preferencial é usada da seguinte forma:

{\ displaystyle G = \ int _ {a} ^ {b} g (x) \, {\ mbox {d}} x = \ int _ {a} ^ {b} {\ frac {g (x)} { f ^ {*} (x)}} f ^ {*} (x) \, {\ mbox {d}} x \ quad = \ quad \ int _ {F (a)} ^ {F (b)} { \ frac {g (F ^ {- 1} (y))} {f ^ {*} (F ^ {- 1} (y))}} \, {\ mbox {d}} y}

onde fizemos a mudança da variável y = F ( x ) com . Parece claro que se então a função a ser integrada à direita está próxima da constante 1 (portanto, de baixa variância). ${\ displaystyle F \, '(x) = f ^ {*} (x)> 0}$ ${\ displaystyle f ^ {*} \ simeq g}$

Para fazer a ligação com a interpretação probabilística da seção anterior, notamos que se define até um fator K que desaparece no quociente. Podemos, portanto, impor isso , o que torna uma densidade de probabilidade em [ a , b ]. A mudança de variável é então interpretada naturalmente como uma mudança de probabilidade e temos a simplificação: $f ^ {*}$ ${\ displaystyle \ int _ {a} ^ {b} f ^ {*} (x) \, {\ mbox {d}} x = 1}$

{\ displaystyle \ int _ {F (a)} ^ {F (b)} {\ frac {g (F ^ {- 1} (y))} {f ^ {*} (F ^ {- 1} ( y))}} \, {\ mbox {d}} y \ quad = \ quad \ int _ {0} ^ {1} {\ frac {g (F ^ {- 1} (y))} {f ^ {*} (F ^ {- 1} (y))}} \, {\ mbox {d}} y}

Essa técnica é imediatamente generalizada em qualquer dimensão.

Aplicação: estimativa de uma probabilidade

Considere que queremos estimar por simulação a probabilidade p t de um evento onde X é uma variável aleatória de distribuição e função de densidade . Esse problema se resume à apresentação geral no sentido de que implementa uma integral a ser estimada. Uma amostra distribuída de forma idêntica e independente (iid) é desenhada nesta lei. Denotamos por k t o número de realizações maior que t . A variável k t é uma variável aleatória seguindo uma distribuição binomial dos parâmetros K e p t : ${\ displaystyle {X \ geq t \}}$ $F$ ${\ displaystyle f (x) = F '(x) \,}$ ${\ displaystyle \ left (X_ {i} \ right) _ {i \ in \ {1, \ dots, K \}}}$

{\ displaystyle P (k_ {t} = k) = {K \ escolha k} p_ {t} ^ {k} (1-p_ {t}) ^ {Kk}, \, \ quad \ quad k = 0, 1, \ pontos, K}

o que significa, em particular, que : a frequência empírica, portanto, converge para sua probabilidade associada p t . ${\ displaystyle \ mathbb {E} (k_ {t}) = Kp_ {t}}$ ${\ displaystyle k_ {t} / K}$

A amostragem preferencial entra em jogo aqui para diminuir a variância da estimativa de Monte-Carlo da probabilidade p t . Na verdade, p t é dado por

{\ displaystyle {\ begin {alinhados} p_ {t} & = {\ mathbb {E}} [I_ {X \ geq t}] \\ & = \ int (I_ {x \ geq t}) {\ frac { f (x)} {f ^ {\ ast} (x)}} f ^ {\ ast} (x) \, dx \\ & = {\ mathbb {E} _ {*}} [(I_ {X \ geq t}) w (X)] \ end {alinhado}}}

onde, nós posamos de novo

{\ displaystyle w (\ cdot) \ equiv {\ frac {f (\ cdot)} {f ^ {\ ast} (\ cdot)}}}

A última igualdade da equação anterior sugere o estimador do seguinte: $p_ {t}$

{\ displaystyle {\ hat {p}} _ {t} = {\ frac {1} {K}} \, \ sum _ {i = 1} ^ {K} (I_ {X_ {i} \ geq t} ) w (X_ {i}), \, \ quad \ quad X_ {i} \ sim f ^ {\ ast}}

É um estimador EP de p t que não tem vieses. Esta sendo definida, o procedimento de estimativa é para gerar uma amostra iid a partir da densidade e para cada realização superior t para calcular o peso W . O resultado será a média obtida com K impressões. A variação deste estimador é: ${\ displaystyle f ^ {\ ast} \,}$

{\ displaystyle {\ begin {align} {\ mbox {Var}} ^ {\ ast} {\ hat {p}} _ {t} & = {\ frac {1} {K}} {\ mbox {Var} } ^ {\ ast} [(I_ {X \ geq t}) w (X)] \\ & = {\ frac {1} {K}} {\ mbox {Var}} ^ {\ ast} [(I_ {X \ geq t}) w (X)] \\ & = {\ frac {1} {K}} \ left [\ mathbb {E} _ {*} [(I_ {X \ geq t}) ^ { 2} w ^ {2} (X)] - p_ {t} ^ {2} \ right] \\ & = {\ frac {1} {K}} \ left [\ mathbb {E} [(I_ {X \ geq t}) ^ {2} w (X)] - p_ {t} ^ {2} \ direita] \ end {alinhado}}}

Aqui, novamente, a densidade deve ser perfilada da melhor forma possível , a fim de reduzir a variância. ${\ displaystyle f ^ {\ ast}}$

Exemplos numéricos

Integração da função beta usando uma distribuição triangular

Detalhe do método

Queremos estimar a seguinte quantidade:

{\ displaystyle G = \ int _ {0} ^ {1} x ^ {4} (1-x) ^ {2} \, {\ mbox {d}} x}

que passa a ser a função beta do parâmetro (5; 3), que é G = 1/105 = 0,0095238095238095. Isso corresponde ao caso geral com a = 0, b = 1 e . ${\ displaystyle g (x) = x ^ {4} (1-x) ^ {2}}$

Simulamos uma amostra de acordo com a lei uniforme padrão para obter o estimador Monte-Carlo clássico: ${\ displaystyle (y_ {1}, \ cdots, y_ {n})}$ ${\ displaystyle {\ mathcal {U}} [0; 1]}$

{\ displaystyle {\ hat {g}} _ {1} = {\ frac {1} {n}} \ sum _ {i = 1} ^ {n} g (y_ {i})}

e o estimador de sua variância:

{\ displaystyle {\ hat {\ sigma}} _ {g_ {1}} ^ {2} = {\ frac {1} {n}} \ left [{\ frac {1} {n}} \ sum _ { i = 1} ^ {n} g ^ {2} (y_ {i}) - {\ hat {g}} _ {1} ^ {2} \ right]}

Inspirando-nos na forma geral da função beta, podemos substituir a distribuição uniforme padrão pela distribuição triangular . ${\ displaystyle {\ mathcal {T}} (a = 0, c = 2/3, b = 1)}$

Parece um triângulo baseado no segmento [0; 1] e "culminando" em (2/3; 2). Sua densidade é

{\ displaystyle f ^ {\ ast} (x) = {\ begin {cases} 3x & x \ in [0; 2/3] \\ 6 (1-x) & x \ in [2/3; 1] \ end {casos}}}

Simulamos uma amostra nesta lei, pelo método da transformada inversa , e, ao definir , o estimador EP é dado por ${\ displaystyle (z_ {1}, z_ {2}, \ cdots, z_ {n})}$ ${\ displaystyle w (x) = g (x) / f ^ {\ ast} (x)}$

{\ displaystyle {\ hat {g}} _ {2} = {\ frac {1} {n}} \ sum _ {i = 1} ^ {n} w (z_ {i})}

e o estimador de sua variância é

{\ displaystyle {\ hat {\ sigma}} _ {g_ {2}} ^ {2} = {\ frac {1} {n}} \ left [{\ frac {1} {n}} \ sum _ { i = 1} ^ {n} w ^ {2} (z_ {i}) - {\ hat {g}} _ {2} ^ {2} \ right]}

Na tabela, vemos que o uso do PE sistematicamente permite reduzir a variância da estimativa em relação à estimativa de Monte-Carlo de mesmo tamanho (ou seja, a n dado). Notamos também que a variância da estimativa é proporcional a 1 / n : indo de n = 1000 para n = 10.000 (multiplicação por 10 do tamanho), reduzimos a variância por um fator de 10.

Comparação do método de Monte-Carlo e amostragem preferencial

não	estimador	tendência	variância	estimador	tendência	variância
	Monte-Carlo Clássico			Amostragem preferencial
500	0,009843	-3,19E-004	1.32E-007	0,009712	-1,88E-004	2.50E-008
1000	0,009735	-2,12E-004	6,53E-008	0,009680	-1,57E-004	1.26E-008
2500	0,009628	-1,04E-004	2.60E-008	0,009576	-5,18E-005	5,36E-009
5000	0,009717	-1,93E-004	1,31E-008	0,009542	-1,83E-005	2.71E-009
10.000	0,009634	-1,10E-004	6,52E-009	0,009544	-1,99E-005	1.35E-009

Espera-se melhorar ainda mais o desempenho considerando uma densidade “mais próxima” da densidade f . O principal problema será fazer simulações. Nos casos mais simples, como a lei triangular, o método da transformação inversa pode ser suficiente; em casos mais complexos, o método de rejeição deve ser usado . ${\ displaystyle f ^ {\ ast}}$

Integração de um Gaussiano

A dificuldade de integrar tal função é que a variável de integração assume valores . Nesse caso, usar uma densidade de probabilidade uniforme é impossível, porque a probabilidade de ocorrência de seria . No entanto, em uma integração Monte-Carlo bruta , não há conhecimento a priori da forma da função, e todos os valores de x devem ter uma probabilidade igual. Assim, a amostragem preferencial é um método que permite integrar funções para uma variável de integração entre a , quando a própria distribuição utilizada permite fornecer valores em . $\ mathbb {R}$ $x \ in {\ mathbb {R}}$ ${\ displaystyle p (x) = \ lim _ {l \ to \ infty} {\ frac {x} {l}} = 0}$ $- \ infty$ $+ \ infty$ $\ mathbb {R}$

Veja também

Links internos

Método Monte-Carlo
(en) J. Morio e M. Balesdent , Estimation of Rare Event Probabilities in Complex Aerospace and Other Systems: uma abordagem prática , Cambridge, Elsevier Science,2015, 216 p. ( ISBN 978-0-08-100091-5 )