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

Regressão Linear Bayesiana com varias variaveis

0 votos
13 visitas
perguntada Jul 6 em Aprendizagem de Máquinas por Victor Candido (21 pontos)  

Resolveremos o seguinte problema do livro Baysian Econometrics de Gary Koops:

Problema 8 PG 57

(a) "Generate an artificial data set of size N D 100 for the Normal linear
regression model with an intercept and one other explanatory variable.
Set the intercept to 0, the slope coefficient to 1:0 and h D 1:0. Generate
the explanatory variable by taking random draws from the U.0; 1/
distribution"

(b) Calculate the posterior mean and standard deviation for this data set using
a Normal-Gamma prior with \(beta=(0,1),v=I_{2,},s^{-2}=1,v=1\)

Compartilhe

1 Resposta

0 votos
respondida Jul 6 por Victor Candido (21 pontos)  
editado Jul 6 por Victor Candido
#Primeiro vamos importar as bibliotecas necessárias
import numpy as np
import pandas as pd
import scipy as sp
import statsmodels as sm
import matplotlib.pyplot as plt
import seaborn as sn

np.random.seed(12345) # Estamos especificando um seed, para que caso alguém queira repetir o exercício, chegue ao mesmo resultado

#Parâmetros iniciais da questão
N = 100
beta1 = 0
beta2 = 1 
h = 1 
sigma = np.sqrt(1/h)

# Vetores Iniciais
x = np.ones((N,2), dtype=np.float64)
y = np.zeros((N,1), dtype=np.float64)

# loop que gera os dados
for i in range(N):
x[i,1] = np.random.uniform(low=0,high=1)
e = sigma* np.random.normal(0,1)
y[i,0] = beta1 + beta2*x[i,1] + e
data = np.concatenate((y,x), axis = 1)

Agora nos vamos calcular a média a posteriori e o desvio padrão para o coeficiente angular(beta2) para esse dataset. Nós faremos isso usando uma distribuição gammar normal com as seguintes especificações: \(beta=(0,1),v=I_{2,},s^{-2}=1,v=1\)

k = x.ndim # número de dimensões para x

Agora vamos inicializar os hiperparâmetros

v_0 = 1
s2inv_0 = 1
s02 = 1/s2inv_0
b0 = np.zeros((k,1), dtype=np.float64)
b0[k-1,0] = 1
c = 1
capv0 = c*np.identity(k, dtype=float)

xsquare = np.dot(x.transpose(), x)
xsquare_inv = np.linalg.inv(xsquare)

Explicitando alguns resultados

bols = (xsquare_inv @ x.transpose()) @ y
s2 = np.transpose(y - x@bols)@(y-x@bols)/(N-k)
v= N - k

Onde \(v=98 e s2= 1,10855\)

Também disponibilizo meu código para quem tiver maior interesse nas respostas e em como fazer: https://drive.google.com/file/d/1rmmSoNB2-5d5GqOokyXj5hN4ZW_vzuPk/view?usp=sharing

...