comparison_of_aic.py

from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression, Lasso, Ridge
from sklearn.pipeline import Pipeline
import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return 2 + 3*x - x**2

X = (np.random.rand(50) * 10 - 5)
Y = f(X) + np.random.normal(0, 5, size=X.shape)
plt.plot(X, Y, 'x')

plt.plot(np.arange(-5, 5, 0.01), f(np.arange(-5, 5, 0.01)), label='true')

get_model = lambda deg:Pipeline([
    ('poly', PolynomialFeatures(degree=deg)),
    ('regression', LinearRegression(fit_intercept=False))])

def plot_model(model, name):
    X = np.arange(-5, 5, 0.01)
    eY = model.predict(X[:, np.newaxis])
    plt.plot(X, eY, '-', label=name)

def likelihood(model, X, y):
    ey = model.predict(X[:, np.newaxis])
    rss = np.linalg.norm(y - ey)**2
    n = y.size
    return -n/2 * (np.log(2*np.pi) + 1 + np.log(rss/n))

def proc(name, model, k):
    model = model.fit(X[:, np.newaxis], Y)

    logL = likelihood(model, X, Y)
    aic = -2 * logL + 2*k
    print('%s model log-likelihood: %8f AIC: %8f' % (name, logL, aic))

    plot_model(model, '%s (AIC=%.1f)' % (name, aic))

for deg in range(5):
    proc('degree=%d' % deg, get_model(deg), deg+1)

plt.ylim(-40, 15)
plt.xlim(-5, 5)
plt.legend(loc='lower right')
plt.savefig('plot.png')
plt.show()

plot.png