scikit-learnで線形回帰

単純な線形回帰

まず簡単な例から。
y = 3 x + 1 + err
という特性を満たすデータからサンプリングを行い、xとyの関係を求める。
ただし、err ~ N(0, 0.1²)とする。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

# create samples
sample_size = 30
err_sigma = 0.1

x = np.random.rand(sample_size, 1)
err = err_sigma*np.random.randn(sample_size, 1)
y = 3 * x + 1 + err

# train a linear regression model
regr = LinearRegression()
regr.fit(x, y)

# make predictions
xt = np.linspace(0.0, 1.0, num=1000).reshape((1000, 1))
yt = regr.predict(xt)

# plot samples and regression result
plt.plot(x, y, 'o')
plt.plot(xt, yt)
plt.show()

二次関数の回帰

次に、多項式基底関数を作って二次関数の回帰を試してみる。

以下の関係を満たすデータからサンプルを行うこととする。
y = -0.3 x² + 1 + err,
err ~ N(0, 0.5²).

scikit-learnにはPolynomialFeaturesという多項式基底を作成する関数があるので、この関数を使えばよい。
また、Pipelineを使うことで基底の作成〜モデルの学習までの処理を纏めることが出来る。

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

# create samples
sample_size = 30
err_sigma = 0.5

x = 10*np.random.rand(sample_size, 1)
err = err_sigma*np.random.randn(sample_size, 1)
y = -0.3 * x*x + 1 + err

# train a linear regression model
regr = Pipeline([
    ('poly', PolynomialFeatures(degree=2)),
    ('linear', LinearRegression())
])
regr.fit(x, y)

# make predictions
xt = np.linspace(0.0, 10.0, num=1000).reshape((1000, 1))
yt = regr.predict(xt)

# plot samples and regression result
plt.plot(x, y, 'o')
plt.plot(xt, yt)
plt.show()

sinの回帰

最後に以下の関係を満たすデータの回帰を行う。

y = sin(x) + err,
err ~ N(0, 0.1²).

基底として用いる多項式の次数を上げることで、よりよくフィットしていく様子が分かる。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
from math import sin

# create samples
sample_size = 100
err_sigma = 0.1

x = 12*np.random.rand(sample_size, 1)
err = err_sigma*np.random.randn(sample_size, 1)
func = np.vectorize(sin)
y = func(x) + err

# plot train data
plt.plot(x, y, 'o')

# train linear regression models with different polynomial basis
deg = [1,3,6]
for d in deg:
    regr = Pipeline([
        ('poly', PolynomialFeatures(degree=d)),
        ('linear', LinearRegression())
    ])
    regr.fit(x, y)

    # make predictions
    xt = np.linspace(0.0, 12.0, num=1000).reshape((1000, 1))
    yt = regr.predict(xt)

    # plot regression result
    plt.plot(xt, yt, label='polynomial of degree %d' % (d))

plt.legend()
plt.show()