import itertools  # Generate combinations
from sklearn.linear_model import LinearRegression
from numpy.random import randn
from scipy import stats
import scipy
import matplotlib.pyplot as plt
import numpy as np

res = LinearRegression()

def RSS_min(X, y, T):
    S_min = np.inf
    m = len(T)
    for j in range(m):
        q = T[j]
        res.fit(X[:, q], y)
        y_hat = res.predict(X[:, q])
        S = np.linalg.norm(y_hat-y)**2
        if S < S_min:
            S_min = S
            set_q = q
    return (S_min, set_q)

df = np.loadtxt("boston.txt", delimiter="\t")
X = df[:, [0, 2, 4, 5, 6, 7, 9, 10, 11, 12]]
p = X.shape[1]
y = df[:, 13]

n, p = X.shape
AIC_min = np.inf
for k in range(1, p+1, 1):
    T = list(itertools.combinations(range(p), k))
    # Generate combinations of k elements from p elements
    S_min, set_q = RSS_min(X, y, T)
    AIC = n*np.log(S_min/n)+2*k
    if AIC < AIC_min:
        AIC_min = AIC
        set_min = set_q
print(AIC_min, set_min)

1619.787608556615 (0, 2, 3, 5, 7, 8, 9)

def IC(X, y, k):
    n, p = X.shape
    T = list(itertools.combinations(range(p), k))
    S, set_q = RSS_min(X, y, T)
    AIC = n*np.log(S/n)+2*k
    BIC = n*np.log(S/n)+k*np.log(n)
    return {'AIC': AIC, 'BIC': BIC}

AIC_seq = []
BIC_seq = []
for k in range(1, p+1, 1):
    AIC_seq.append(IC(X, y, k)['AIC'])
    BIC_seq.append(IC(X, y, k)['BIC'])
x_seq = np.arange(1, p+1, 1)
plt.plot(x_seq, AIC_seq, c="red", label="AIC")
plt.plot(x_seq, BIC_seq, c="blue", label="BIC")
plt.xlabel("Number of variables")
plt.ylabel("AIC/BIC value")
plt.title("Change in AIC and BIC values with the number of variables")
plt.legend()

<matplotlib.legend.Legend at 0x7f6ffae22ef0>

Chapter 5 Information Criteria¶

46¶

48¶