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