Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como utilizar o Python para ilustrar a aplicação da teoria Bayesiana sobre o caso mais simples de regressão linear?

+1 voto
56 visitas

1 Resposta

+1 voto
respondida Jul 1 por Pedro Antero (26 pontos)  
editado Jul 6 por Pedro Antero

Uma forma interessante de respondermos a essa pergunta é combinarmos a teoria apresentada no Capítulo 2 do livro Bayesian Econometrics, de Gary Koop, com uma demonstração computacional em Python.

INTRODUÇÃO TEÓRICA

O modelo mais simples de regressão linear (sem intercepto e com apenas uma variável explicativa) pode ser descrito pela equação abaixo, para \(i=1...N\) observações:

\( y_i=\beta x_i + \epsilon_i \tag{1}\)

Para ilustrar como podemos aplicar a teoria Bayesiana sobre modelos desse tipo, vamos adotar mais 2 hipóteses simplificadoras:

  1. \(\epsilon_i\) é i.i.d. \(N(0, \sigma^2)\)
  2. \(x_i\) é independente de \(\epsilon_i\), para \(i=1...N\)

Para tornar a aplicação da teoria Bayesiana ainda mais simples, utilizaremos o conceito de priors naturais conjugadas. Sabemos que a função da prior é agregar ao modelo quaisquer informações que o pesquisador já tenha de antemão, antes mesmo de observar os dados. Portanto, as priors podem ter qualquer formato. Entretanto, as chamadas priors naturais conjugadas possuem duas vantagens que justificam sua utilização:

  1. Possuem a mesma forma funcional que a função de verossimilhaça; e
  2. Quando combinadas com a função de verossimilhança, geram uma posterior pertencente a mesma classe de distribuição

Usando as hipóteses simplificadoras para \(\epsilon_i\) e \(x_i\), acima, e denotando \(h=\frac{1}{\sigma^2}\), é possível escrever a função de verossimilhança do modelo de regressão linear da Equação (1), como:

\(p(y | \beta, h) = \frac{h^{1/2}}{(2 \pi)^{N/2}}exp\left[-\frac{h}{2}(\beta - \hat{\beta})^2 \sum^N_{i=1}x_i^2\right]\left[h^{\nu/2}exp\left(-\frac{h\nu}{2s^{-2}}\right)\right] \tag{2} \)

Onde \(\hat{\beta}\), \(s\), e \(\nu\) são, respectivamente, os estimadores de \(\beta\) e do desvio padrão (por Mínimos Quadrados Ordinários) e os graus de liberdade tipicamente utilizados na abordagem frequentista:

\(\hat{\beta} = \frac{\sum x_i y_i}{\sum x_i^2} \tag{3} \)

\(s^2 = \frac{\sum({y_i-\hat{\beta}x_i})^2}{\nu} \tag{4} \)

\(\nu = N -1 \tag{5}\)

A dedução da Equação (2) é resultado direto da aplicação das hipóteses sobre os erros, \(\epsilon_i\), e sobre as variáveis explicativas, \(x_i\), e do emprego da densidade de probabilidade da função Normal. Maiores detalhes podem ser encontrados no Capítulo 2 e no Apêndice B, do livro Bayesian Econometrics, de Gary Koop.

Consideremos agora as duas características de uma prior natural conjugada, discutidas acima. É fácil perceber que o formato da Equação (2), sugere que a prior natural conjugada para o modelo da Equação (1) deve envolver necessariamente uma distribuição Normal para \(\beta|h\) e uma distribuição Gama para \(h\). De fato, é possível mostrar que \(\beta|h \sim N(\underline{\beta}, h^{-1}\underline{V}) \) e que \(h \sim G(\underline{s}^{-2}, \underline{\nu}) \). Assim, a prior conjugada natural pode ser denotada pela seguinte função Gama Normal:

\(\beta, h \sim NG(\underline{\beta}, \underline{V}, \underline{s}^{-2}, \underline{\nu}) \tag{6} \)

Onde \((\underline{\beta}, \underline{V}, \underline{s}^{-2}, \underline{\nu})\) são os chamados hiper-parâmetros da prior e refletem as informações que o pesquisador já saiba de antemão sobre o processo em análise.

Sigamos agora ao estudo da posterior. Sabemos que a densidade da posterior é resultado da combinação entre as informações já sabidas pelo pesquisador sobre o processo antes da observação dos dados (prior) e àquelas efetivamente extraídas das observações. Assim, a posterior será proporcional ao produto entre a prior (Equação (6)) e a verossimilhança, dada por \( p(y_i | \beta, \sigma^2) = \frac{1}{(2 \pi)^{N/2}\sigma^N}exp\left[ -\frac{1}{\sigma^2} \sum^N_{i=1} (y_i - \beta x_i)^2 \right] \).

É fácil perceber que essa combinação também resultará em uma função Gama Normal. Assim, a distribuição da posterior terá a forma:
\(\beta, h | y \sim NG (\bar{\beta}, \bar{V}, \bar{s}^{-2}, \bar{\nu}) \tag{7} \)

Onde:

\( \bar{V} = \frac{1}{\underline{V}^{-1} + \sum x_i^2} \tag{8} \)

\( \bar{\beta} = \bar{V}(\underline{V}^{-1} \underline{\beta} + \hat{\beta} \sum x_i^2) \tag{9} \)

\( \bar{\nu} = \underline{\nu} + N \tag{10} \)

e \( \bar{s}^{-2} \) é definido implicitamente por:

\( \bar{\nu s}^2 = \underline{\nu s^2} + \nu s^2 + \frac{(\hat{\beta} - \underline{\beta})^2}{\underline{V} + \frac{1}{\sum x_i^2}} \tag{11}\)

A partir da Equação (7), não é difícil mostrar que as distribuições marginais, para \( \beta | y \) e \( h | y\) serão dadas por:

\( \beta | y \sim t(\bar{\beta}, \bar{s}^2\bar{V}, \bar{\nu}) \tag{12} \) e
\( h| y \sim G(\bar{s}^{-2}, \bar{\nu}) \tag{13} \)

De (11), segue que:

\( E(\beta | y) = \bar{\beta} \tag{14} \)
e
\( var(\beta|y) = \frac{\bar{\nu s}^2}{\nu - 2}\bar{V} \tag{15} \)

De forma análoga, a partir de (12), temos que:

\( E(h | y) = \bar{s}^{-2} \tag{16} \)
e
\( var(h|y) = \frac{2\bar{s}^{-2}}{\bar{\nu}} \tag{17} \)

Abaixo, discutiremos dois tópicos complementares: a comparação entre modelos (MC) e predição (P).

Comparação entre modelos (MC)

Agora, vamos conversar um pouco sobre estratégias de comparação entre duas (ou mais) especificações concorrentes para o processo em análise.

Vamos supor, por exemplo, que temos dois modelos concorrentes, \( M_1 \) e \( M_2 \) que se propõem a explicar \( y \). Ou seja, tomemos \( M_j \), para \( j = 1, 2 \), de modo que, para cada um dos modelos, possamos escrever a Equação (1), como:

\( y_i = \beta_j x_{ji} + \epsilon_{ji} \tag{MC-1} \)

Suponha, ainda, que as hipóteses sobre a independência entre os erros e entre os erros e as variáveis explicativas continuam válidas.

Tendo em mente que a ideia de que a posterior é formada por uma combinação entre a prior e os dados, podemos definir a seguinte razão de probabilidade:
\( PO_{12} = \frac{p(y | M_1)p(M_1)}{p(y | M_2)p(M_2)} \tag{MC-2} \)

Mantendo-se as hipóteses de independência e considerando que a prior a ser utilizada será do tipo natural conjugada, é possível mostrar que:

\( p(y|M_j) = c_j \left( \frac{\bar{V}_j}{\underline{V}_j} \right)^{\frac{1}{2}} (\bar{\nu}_j \bar{s}_j^2)^{-{\frac{\bar{\nu}_j}{2}}} \tag{MC-3} \)

Com \( j = 1,2 \), onde:
\( c_j = \frac{\Gamma \left( \frac{\bar{\nu}_j}{2} \right) (\underline{\nu}_j \underline{s}_j^2)^{-\frac{\underline{\nu}_j}{2}} }{\Gamma \left( \frac{\underline{\nu}_j}{2} \right) \pi^{\frac{N}{2}}} \tag{MC-4} \)

Portanto, para o caso de dois modelos, \(M_1\) e \(M_2\), podemos reescrever a Equação (MC-2) como:

\( PO_{12} = \frac{c_1 \left( \frac{\bar{V}_1}{\underline{V}_1} \right)^{\frac{1}{2}} (\bar{\nu}_1 \bar{s}_1^2)^{-{\frac{\bar{\nu}_1}{2}}} p(M_1)}{c_1 \left( \frac{\bar{V}_2}{\underline{V}_2} \right)^{\frac{1}{2}} (\bar{\nu}_2 \bar{s}_2^2)^{-{\frac{\bar{\nu}_2}{2}}}p(M_2)} \tag{MC-5} \)

Sabendo ainda que \(p(M_j|y)\) devem somar 1 (para \( j = 1,2\)), podemos escrever:
\( p(M_1|y) = \frac{PO_{12}}{1+PO_{12}} \tag{MC-6} \)
e;
\( p(M_2|y) = \frac{1}{1+PO_{12}} \tag{MC-7} \)

A Equação (MC-5) fornece 4 intuições sobre a comparação entre \(M_1\) e \(M_2\)

  1. Quanto maior a razão de probabilidade entre as priors, \( \frac{p(M_1)}{p(M_2)}\), maior será a dominância do modelo \( M_1\);
  2. Sabendo que \( \bar{\nu}_j\bar{s}^2_j \) contém o termo \( \nu_j s^2_j \) (vide Equação (11) ), deduz-se que a razão de probabilidade dará um maior peso ao modelo com menor soma dos quadrados dos erros;
  3. Sabendo que \( \bar{\nu}_j\bar{s}^2_j \) contém o termo \((\hat{\beta}_j - \underline{\beta}_j)^2\) (vide, novamente a Equação (11)), pode-se deduzir, também, que a razão de probabilidade atribuirá um peso maior ao modelo que apresentar a maior coerência entre o estimador da prior e dos dados; e
  4. Tomando o termo \( \frac{\bar{V}_1}{\underline{V}_1} \), conclui-se que a razão de probabilidade recompensará priors com menores variâncias.

Predição (P)

Nesta seção, vamos discutir como utilizar a abordagem Bayesiana para realizar previsões sobre um \( y^*\) não observado. Ou seja, tomemos o modelo abaixo, em que observamos apenas \(x^*\):

\( y_i=\beta x_i + \epsilon_i \tag{P-1}\)

Mantendo as hipóteses de independência dos erros, é possível mostrar que \( y^{*} | y \) segue uma distribuição t de Student da forma:
\( y^{*} | y \sim t(\bar{\beta x^{*}}, \bar{s}^2\{1 + \bar{V}x^{*2}\}, \bar{\nu}) \tag{P-2}\)

IMPLEMENTAÇÃO EM PYTHON

Agora, seguindo a seção 2.7 do livro Bayesian Econometrics, de Gary Koop, vamos tentar implementar computacionalmente alguns dos conceitos estudados na Introdução Teórica.

Primeiro, vamos definir o modelo da Equação (1) com \( N = 50\); \( \beta = 1,5 \); \( x_i\) i.i.d \( N(0, 0,5) \) e \( \epsilon_i\) i.i.d. \( N(0,1) \), ou seja, com \(h=1\).

import numpy as np  # imports numpy library
N = 50  # number of observations
beta = 2
x = np.random.normal(0, 1, N)  # x_i's
epsilon = np.random.normal(0, 1, N)  # epsilon_i's

y = beta*x + epsilon  # Eq. (1)

Consideremos ainda a análise com 2 tipos de priors:

  1. Prior não-informativa
  2. Prior informativa (com \( \underline{\beta} = 1,5; \underline{V} = 0,25; \underline{\nu} = 10\) e \(\underline{s}^{-2}=1\))

Caso 1: Prior não informativa

Uma prior é dita não informativa quando não adiciona nenhuma informação a posterior. Esse tipo de prior é utilizado quando não temos nenhuma pista da distribuição das variáveis antes da observação dos dados.

A partir das Equações (8)-(11) é fácil perceber que para eliminarmos a influência da prior sobre os parâmetros da posterior, basta fazermos \( \underline{\nu} = 0 \) e \( \underline{V} \rightarrow \infty \).

Desta forma, podemos reescrever as Equações (8)-(11) como:

\( \bar{V} = \frac{1}{\sum x_i ^2} \tag{18} \)

\( \bar{\beta} = \hat{\beta} \tag{19} \)

\( \bar{\nu} = N \tag{20} \)

\( \bar{\nu s}^2 = \nu s^2 \tag{21} \)

Observe que, nesse caso, os parâmetros da posterior, acima, dependem apenas da informação observada nos dados. Ou seja, para o caso da prior não informativa, os parâmetros estimados são os mesmos daqueles que seriam obtidos por uma abordagem frequentista em Mínimos Quadrados Ordinários.

A partir das Equações (3), (4), (5), (14), (15), (16), (17), (19) e (21) podemos, então, computar a média e o desvio padrão dos parâmetros da posterior não informativa \( \bar{\beta} \) e \( \bar{h} \).

nu = N - 1   # Eq. (5)
beta_hat = sum(x * y) / sum(x * x)   # Eq. (3)
s_sqd = sum((y - (beta_hat*x)) * (y - (beta_hat*x))) / nu   # Eq. (4)

V_bar_nip = 1 / sum(x * x)
beta_bar_nip = beta_hat  # Eq. (19)
nu_bar_nip = N  # Eq. (20)
s_bar_sqd_nip = (nu * s_sqd) / nu_bar_nip  # Eq. (21)

e_beta_post_nip = beta_bar_nip  # Eq. (14)
var_beta_post_nip = ((nu_bar_nip * s_bar_sqd_nip) / (nu - 2)) * V_bar_nip # Eq. (15)
std_dev_beta_post_nip = np.sqrt(var_beta_post_nip)

e_h_post_nip = 1 / s_bar_sqd_nip # Eq. (16)
var_h_post_nip = (2 * (1 / s_bar_sqd_nip)) / nu_bar_nip # Eq. (17)
std_dev_h_post_nip = np.sqrt(var_h_post_nip)

print('e_beta_post_nip')
print(e_beta_post_nip)
print('std_dev_beta_post_nip')
print(std_dev_beta_post_nip)
print('e_h_post_nip')
print(e_h_post_nip)
print('std_dev_h_post_nip')
print(std_dev_h_post_nip)

Caso 2: Prior informativa

A média e o desvio padrão dos parâmetros da prior informativa podem ser calculados por meio das equações (14)-(17). Devemos, entretanto, ter a atenção de empregarmos os parâmetros referentes a prior, já mencionados acima: \( \underline{\beta} = 1,5; \underline{V} = 0,25; \underline{\nu} = 10\) e \(\underline{s}^{-2}=1\).

# Hyper-parameters of the informative prior
beta_uline = 1.5
V_uline = 0.25
nu_uline = 10
s_minus_sqd_uline = 1

e_beta_prior = beta_uline  # Eq. (14)
var_beta_prior = ((nu_uline * (1 / s_minus_sqd_uline)) / (nu_uline - 2)) * V_uline  # Eq. (15)
std_dev_beta_prior = np.sqrt(var_beta_prior)

e_h_prior = s_minus_sqd_uline  # Eq. (16)
var_h_prior = (2 * s_minus_sqd_uline) / nu_uline  # Eq. (17)
std_dev_h_prior = np.sqrt(var_h_prior)

print('e_beta_prior')
print(e_beta_prior)
print('std_dev_beta_prior')
print(std_dev_beta_prior)
print('e_h_prior')
print(e_h_prior)
print('std_dev_h_prior')
print(std_dev_h_prior)

Para estimarmos a média e o desvio padrão dos parâmetros \( \bar{\beta} \) e \( \bar{h} \) da posterior informativa, basta utilizarmos as Equações (3)-(5), (8)-(11) e (14)-(17) e os parâmetros \( \underline{\beta}, \underline{V}, \underline{\nu}\) e \(\underline{s}^{-2}\), já empregados acima.

V_bar_ip = 1 / ((1 / V_uline) + sum(x * x))  # Eq. (8)
beta_bar_ip = V_bar_ip * (((1 / V_uline) * beta_uline) + (beta_hat * sum(x * x)))  # Eq. (9)
nu_bar_ip = nu_uline + N  # Eq. (10)
s_bar_sqd_ip = (1 / nu_bar_ip) * ((nu_uline * (1 / s_minus_sqd_uline)) + (nu * s_sqd) + ((beta_hat - beta_uline) ** 2 / (V_uline + (1/(sum(x * x))) )))  # Eq. (11)

e_beta_post_ip = beta_bar_ip  # Eq. (14)
var_beta_post_ip = ((nu_bar_ip * s_bar_sqd_ip) / (nu_bar_ip - 2)) * V_bar_ip  # Eq. (15)
std_dev_beta_post_ip = np.sqrt(var_beta_post_ip)

e_h_post_ip = 1 / s_bar_sqd_ip  # Eq. (16)
var_h_post_ip = (2 * (1 / s_bar_sqd_ip)) / nu_bar_ip  # Eq. (17)
std_dev_h_post_ip = np.sqrt(var_h_post_ip)

print('e_beta_post_ip')
print(e_beta_post_ip)
print('std_dev_beta_post_ip')
print(std_dev_beta_post_ip)
print('e_h_post_ip')
print(e_h_post_ip)
print('std_dev_h_post_ip')
print(std_dev_h_post_ip)

Os resultados das propriedades dos parâmetros obtidas a partir do código acima são mostrados nas Tabelas 1 e 2 abaixo para \( \beta \) e \( h \), respectivamente.

A imagem será apresentada aqui.
Propriedades do parâmetro \( \beta \)

A imagem será apresentada aqui.
Propriedades do parâmetro \( h \)

O Gráfico 1 representa as priors e posteriors marginais para \( \beta\). Ele pode ser construído a partir do seguinte código:

x = np.linspace(e_beta_prior - 0.5, e_beta_prior + 1.5, 100)

fig, ax = plt.subplots()

line1, = ax.plot(x, t.pdf(x, nu_uline, beta_uline, np.sqrt(var_beta_prior)), 'k--', lw=1, alpha=1, label='Prior')
line2, = ax.plot(x, t.pdf(x, nu_bar_ip, e_beta_post_ip, np.sqrt(var_beta_post_ip)), 'k-', lw=1, alpha=1, label='Posterior')
line3, = ax.plot(x, t.pdf(x, nu_bar_nip, e_beta_post_nip, np.sqrt(var_beta_post_nip)), 'k-.', lw=1, alpha=1, label='Likelihood')

ax.legend(loc='upper right')
plt.savefig('marg_pdf.eps')

plt.show()

A imagem será apresentada aqui.
Gráfico 1 - Priors e posteriors marginais para \( \beta \)

Conforme previsto pela Equação (12), as curvas seguem uma distribuição t de Student. Observe também que a posterior calculada com base na prior não informativa se confunde com a curva de "verossimilhança" (Likelihood).

As tabelas e o gráfico acima ilustram como as informações da prior e dos dados são combinadas para a formação da "posterior". No gráfico, por exemplo, percebe-se claramente que a curva da posterior, em linha sólida, representa um "meio termo" entre as curvas pontilhadas da prior e da verossimilhança. A mesma conclusão pode ser retirada dos parâmetros apresentados nas tabelas. Os resultados se assemelham bastante àqueles da seção 2.7 do livro Bayesian Econometrics, de Gary Koop. As pequenas diferenças observadas se devem à distribuição aleatória das variáveis.

Comparação entre modelos

Para ilustrarmos a comparação entre modelos, iremos testar o modelo com que vínhamos trabalhando até aqui (\(M_1\)) contra uma especificação que contém somente um intercepto, ou seja, para a qual \( x_i = 1 \) para \( i= 1, ..., 50\), que chamaremos de \(M_2\).

Para simplificarmos os cálculos, vamos assumir que ambos os modelos utilizam a mesma prior informativa utilizada anteriormente (\( NG(1,5, 0,25, 1, 10) \)). Ou seja, \( \frac{p(M_1)}{p(M_2)} = 1\)

O código abaixo calcula a razão de probabilidades da posterior utilizando as Equações (MC-4) e (MC-5).

c_post_ip = special.gamma(nu_bar_ip / 2) * ((nu_uline * (1 / s_minus_sqd_uline)) ** (nu_uline / 2)) / (special.gamma(nu_uline / 2) * (np.pi ** (N / 2))) # Eq. (MC-4)
PO12_num = c_post_ip * (( V_bar_ip / V_uline ) ** (1 / 2)) * ((nu_bar_ip * s_bar_sqd_ip) ** (- nu_bar_ip / 2))

x2 = np.ones((N, 1))

beta_hat_2 = np.dot(np.dot(np.linalg.inv(np.dot(x2.transpose(), x2)), x2.transpose()), y)
s_sqd_2 = np.dot(np.transpose(y - np.dot(x2, beta_hat_2)), y - np.dot(x2, beta_hat_2)) / nu

V_bar_2 = 1 / ((1/(V_uline) + np.dot(x2.transpose(), x2)))

nu_bar_2 = nu_bar_ip

s_bar_sqd_2 = (1 / nu_bar_2) * ((nu_uline * (1 / s_minus_sqd_uline)) + (nu * s_sqd_2) + (((beta_hat_2 - beta_uline) ** 2 )/(V_uline + np.linalg.inv(np.dot(x2.transpose(), x2)))))  # Eq. (11)
c_2 = special.gamma(nu_bar_2 / 2) * ((nu_uline * (1 / s_minus_sqd_uline)) ** (nu_uline / 2)) / (special.gamma(nu_uline / 2) * (np.pi ** (N / 2)))
PO12_den = c_2 * (( V_bar_2 / V_uline ) ** (1 / 2)) * ((nu_bar_2 * s_bar_sqd_2) ** (- nu_bar_2 / 2))

PO12 = PO12_num / PO12_den # Eq. (MC-5)

p_M1_y = PO12 / (1 + PO12) # Eq (MC-6)

Conforme esperado, o resultado fornece \( PO_{12} \approx 1 \). Isso indica que o modelo \( M_1\), construído com base no processo real de geração dos dados da Equação (1) tem quase 100% de chances de explicar melhor \( y \) que o modelo \( M_2 \), que inclui apenas um intercepto.

Predição

Para ilustrar a teoria sobre predição, vamos supor que queiramos saber o comportamento do \( y^*\), não observável, associado à \( x^* = 0,5 \).

Usando a prior não informativa, temos que a equação (P-2) fica:
\( y^{*} | y \sim t(1,05, 0,86, 50) \)

Já com a prior informativa, temos que \( y^{*} | y\) seguirá a seguinte distribuição:
\( y^{*} | y \sim t(0,92, 0,96, 60) \)

Os códigos abaixo nos ajudam a encontrar os parâmetros das distribuições t de Student, acima.

x_star = 0.5

print('beta_bar_nip * x_star')
print(beta_bar_nip * x_star)

print('(s_bar_sqd_nip 2) * (1 + V_bar_nip * (x_star ** 2))')
print((s_bar_sqd_nip) * (1 + V_bar_nip * (x_star ** 2)))

print('nu_bar_nip')
print(nu_bar_nip)

print('beta_bar_ip * x_star')
print(beta_bar_ip * x_star)

print('(s_bar_sqd_ip) * (1 + V_bar_nip * (x_star ** 2))')
print((s_bar_sqd_ip) * (1 + V_bar_ip * (x_star ** 2)))

print('nu_bar_ip')
print(nu_bar_ip)
comentou Jul 7 por IgorNascimento (51 pontos)  
Ficou claro o processo de estimação do modelo de regressão descrito utilizando inferência bayesiana. É importante salientar que o mais interessante da abordagem bayesina é a possibilidade do processo "on-line" de estimação. A medida que se tenha acesso a uma novas informações \(x_i,y_i\), os parâmetros \(\beta\) e \(\sigma\) são reestimados. Isto é comum quando as variáveis de estudo são séries temporais. Nesse caso, a posteriori obtida no passo \(t\) é a nova "priori informativa" para o processo de estimação em \(t+1\).

Além disso, o processo de MC é útil quando o modelo regressão não permite posteriori com forma fechada, isto é, solução analítica.
...