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

Como implementar o algorítmo de seleção de variáveis conhecido como Forwards Stagewise em Python?

+1 voto
72 visitas
perguntada Jul 7, 2016 em Aprendizagem de Máquinas por Matheus Costa (26 pontos)  

Considere a questão 5.18 do livro Modern Multivariate Statistical Techniques - Alan Julian Izenman que versa sobre o uso da base de dados "Boston housing data"e variable selection algorithms. Interprete os resultados e documente cuidadosamente o seu material e implementação. Como implementar o algoritmo de seleção de variáveis em Python conhecido como Forwards Stagewise em Python? Compare esses resultados com as estimativas de OLS usando os dados descritos acima.

Compartilhe

1 Resposta

0 votos
respondida Jul 8, 2016 por Matheus Costa (26 pontos)  
editado Jul 8, 2016 por Matheus Costa

O modelo, como já salientado na questão, faz parte da categoria de algoritmos de escolha de variáveis, no sentido em que inclui as variáveis no modelo desenvolvida na medida em que estas mostram um potencial explicativo dentro da previsão. O Forwards Stagewise, diferentemente, por exemplo, do Forwards Stepwise, se utiliza de um processo iterativo não para a escolha per se das variáveis, mas para a construção dos parâmetros de \( \beta \). Fora isso o modelo funciona como um MQO tradicional. Faremos a implementação, como sugerido, seguindo a seção 5.9 do livro de Alan Izenman citado na pergunta, utilizando, na medida do possível.

O primeiro passo para a implementação é normalizar o vetor de variáveis \(X\), de modo que tenha média zero e variância um, e a variável independente de modo a garantir média zero. Para isso definimos uma função de normalização para \(X\) e em seguida subtraímos \(y\) de sua média:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
boston = load_boston()
            # Dados obtidos da biblioteca sklearn

y= boston.target
X= np.array(boston.data)
i = np.array([[1] for j in range(len(X))])
n = len(X)
k = len(X[1])

print X.shape, n, k

# Normalizando para que o vetor de variaveis
# tenha media zero e norma um.

def Norm(X):
    Xs = np.zeros(k)
    Xn = np.zeros((n,k))
    Xm = np.sum(X,axis = 0)/n
    for i in range(n):
        Xn[i] = X[i] - Xm
    for i in range(k):
        for j in range(n):
            Xs[i] = Xs[i] + Xn[j][i]**2
    for i in range(k):
        for j in range(n):
            Xn[j][i] = Xn[j][i]/np.sqrt(Xs[i])
    return Xn

X = Norm(X)
y = y - sum(y)/n
SST = sum(np.power(y,2))

Em seguida definimos a função para o modelo conforme o passo-a-passo descrito no livro citado, aqui levemente modificado para se adequar melhor à implementação:

1 - Inicialize o vetor \( \beta = 0 \), tal que \( X\beta = 0 \) e \( r = y \).

2 - Encontre a variável de \( X \), digamos \( j_1 \) que possuí maior correlação com \( r \). Ou seja, \( j_1 = argmax_j(|cov(X_j,r)|)\)

3 - Atualize \( \beta_{j_1} \) para \( \beta_{j_1} + \epsilon*sign(cov(X_j,r)) \).

4 - Atualize \( r \) para \( r = y - X\beta \), usando o novo valor de \(\beta\).

5 - Repita os passos 2 - 4 até que \(|cov(X,r)|\) esteja próximo do vetor de zeros.

def Fstage(X,y,eps,steps,color):
    i = 0  #iterador
    r = y[:]
    b = np.zeros(k) #beta
    plot = np.zeros((steps,2)) #vetor para plotar
    while i < steps:
        cov = np.matmul(r,X)
        j1 = np.argmax(np.absolute(cov))
        b[j1] = b[j1] + np.sign(cov[j1])*eps
        r = y - np.matmul(b,X.T)
        SSR = np.dot(r,r) #Soma dos quadrados dos resíduos
        plot[i] = [i,np.absolute(cov[j1])]
        i = i + 1
    plt.plot(plot[:,0],plot[:,1],color,markersize=4)

    SSR = np.dot(r,r)

    return np.array(b), 1-SSR/SST

Com finalidade de comparação, implementamos também um MQO.

def linreg(X,y):
                # Implementacao da regressao linear
    n=len(X)
    Xt = X.T
    XXI=np.linalg.inv(np.matmul(Xt,X))
    Xy = np.matmul(Xt,y)
    b = np.matmul(XXI,Xy)
    res = [0 for i in range(n)]
    SSR = 0
    for i in range (n):
        res[i] = y[i] - np.matmul(X,b)[i]
        SSR = SSR + res[i]**2
    R2 = 1 - SSR/SST
    return np.array(b), R2

Código para plotar os gráficos:

plt.xlabel('Passos', fontsize=14)
plt.ylabel('SSR', fontsize=14)
plt.text(200, 140, r'$\epsilon = 0.5$', color ='y', fontsize=16)
plt.text(200, 130, r'$\epsilon = 1$', color ='b',fontsize=16)
plt.text(200, 120, r'$\epsilon = 5$', color ='r',fontsize=16)

E por fim, imprimindo os resultados:

Comp = np.zeros((k,2))
for i in range(k):
    Comp[i] = [linreg(X,y)[0][i],Fstage(X,y,1,500,u'b-')[0][i]]

print 'Comparacao entre MQO e Fstage com $\epsilon = 1$', Comp
print 'R2 MQO:', linreg(X,y)[1]
print 'R2 para $\epsilon = 5$:',Fstage(X,y,5,500,u'r-')[1]
print 'R2 para $\epsilon = 1$:',Fstage(X,y,1,500,u'b-')
print 'R2 para $\epsilon = 0.5$:', Fstage(X,y,0.5,500,u'y-')[1]
print 'Coeficientes para $\epsilon = 1$', Fstage(X,y,1,500,u'b-')[0]
plt.show()

De modo que obtemos:

Comparacao entre MQO e Fstage com \(\epsilon = 1\) [
[-20.7041 | -20. ]
[ 24.3160 | 23. ]
[ 3.21596 | 1. ]
[ 15.3457 | 16. ]
[-46.3406 | -45. ]
[ 60.0745 | 61. ]
[ 0.47509 | 0. ]
[-69.8328 | -69. ]
[ 59.8079 | 56. ]
[-46.6961 | -42. ]
[-46.3870 | -46. ]
[ 19.2696 | 19. ]
[-84.3244 | -84. ]]
R2 MQO: 0.740607742865
R2 para \(\epsilon = 5\): 0.740383500337
R2 para \(\epsilon = 1\): 0.740506010253
R2 para \(\epsilon = 0.5\): 0.70764037213

A imagem será apresentada aqui.

Fácil perceber, tanto pelos dados quanto pela imagem, que para quaisquer que sejam os valores de \(\epsilon\) a convergência para os coeficientes de MQO ocorre, ainda que seja necessária uma maior quantidade de passos para isso. Além disso, valores elevados de \(\epsilon \) impõe uma restrição à capacidade de adaptação dos coeficientes, como consequência a capacidade de previsão fica restrita e ao se aproximar dos valores de convergência os parâmetros se mantém constantemente variando.

...