t.py

# -*- coding: utf-8 -*-

import numpy as np
from math import log, pi
from sklearn.linear_model import LinearRegression

N = 100
M = 5
NUM_TRIAL = 1000
before = np.zeros(M)
after = np.zeros(M)

def true_model(X):
    first_column = X[:, 0]
    T = first_column + np.random.randn(N)
    return T

def gaussian_loglik(x, mu, sigma=1):
    return (
        - log(2 * pi * sigma ** 2) / 2
        - ((x - mu) ** 2).sum() / (2 * sigma ** 2))

for i in range(NUM_TRIAL):
    # train
    X = np.random.randn(N, M)
    T = true_model(X)
    lrs = [LinearRegression() for j in range(M)]
    for j in range(M):
        data = X[:, :j + 1]  # take j columns from X
        lrs[j].fit(data, T)
        predict = lrs[j].predict(data)
        before[j] += gaussian_loglik(predict, T)

    # test
    X = np.random.randn(N, M)
    T = true_model(X)
    for j in range(M):
        data = X[:, :j + 1]  # take j columns from X
        predict = lrs[j].predict(data)
        after[j] += gaussian_loglik(predict, T)

print (before - before[0]) / NUM_TRIAL
#[ 0.          0.49161953  0.99283812  1.47318384  1.9262704 ]

print (after - after[0]) / NUM_TRIAL
#[ 0.         -0.55409662 -1.11481264 -1.66386383 -2.08548524]