概要

あえて重回帰を Chainer でやってみました。説明変数が2個だけの最小構成の重回帰ですが。

重回帰

この記事からデータをお借りしました。ただ、単位を直しました。

いま、下表のようなデータがあるとします。

身長 $x_1$ (m)	ウェスト $x_2$ (m)	体重 $y$ (100kg)
1.65	0.65	0.50
1.70	0.68	0.60
1.72	0.70	0.65
1.75	0.65	0.65
1.70	0.80	0.70
1.72	0.85	0.75
1.83	0.78	0.80
1.87	0.79	0.85
1.80	0.95	0.90
1.85	0.97	0.95

身長 $x_1$ とウェスト $x_2$ を説明変数として、体重 $y$ を予測します。

コード

multi_regression.py

import numpy as np
import chainer
from chainer import Chain, Variable
import chainer.functions as F
import chainer.links as L

xs = np.array(
    [[1.65, 0.65],
     [1.70, 0.68],
     [1.72, 0.70],
     [1.75, 0.65],
     [1.70, 0.80],
     [1.72, 0.85],
     [1.83, 0.78],
     [1.87, 0.79],
     [1.80, 0.95],
     [1.85, 0.97]]
    , "f")
ys = np.array(
    [[0.50],
     [0.60],
     [0.65],
     [0.65],
     [0.70],
     [0.75],
     [0.80],
     [0.85],
     [0.90],
     [0.95]]
    , "f")

class MultiRegression(Chain):
    def __init__(self):
        super().__init__(
            l1=L.Linear(2, 1)
        )
    def __call__(self, x):
        return self.l1(x)

def normal_equation(xs, ys):
    xs2 = np.hstack((np.ones((xs.shape[0], 1), "f"), xs))
    theta = np.linalg.inv(xs2.T.dot(xs2)).dot(xs2.T).dot(ys)
    return theta

def normalize(xs, ys):
    xmeans = np.mean(xs, axis = 0)
    print("xmeans = " + np.array_str(xmeans))
    xstds  = np.std(xs, axis = 0)
    print("xstds = " + np.array_str(xstds))
    ymean = np.mean(ys)
    print("ymean = %f" % ymean)
    ystd  = np.std(ys)
    print("ystd = %f" % ystd)
    nxs = (xs - xmeans) / xstds
    nys = (ys - ymean) / ystd
    return nxs, nys, xmeans, xstds, ymean, ystd

def denormalize_params(nW, nb, xmeans, xstds, ymean, ystd):
    W = ystd / xstds * nW
    b = ymean - (sum(ystd / xstds * xmeans * nW) + ystd * nb)
    return W, b

def train(xs, ys):
    xs = Variable(xs)
    ys = Variable(ys)

    model = MultiRegression()
    alpha = 0.5
    optimizer = chainer.optimizers.SGD(lr = alpha)
    optimizer.setup(model)

    for i in range(100):
        model.l1.zerograds()
        yp = model(xs)
        loss = F.mean_squared_error(yp, ys)
        loss.backward()
        print("=== Epoch %d ===" % (i + 1))
        print("loss = %f" % loss.data)
        print("model.l1.W.data = %s" % np.array_str(model.l1.W.data))
        print("model.l1.W.grad = %s" % np.array_str(model.l1.W.grad))
        print("model.l1.b.data = %s" % np.array_str(model.l1.b.data))
        print("model.l1.b.grad = %s" % np.array_str(model.l1.b.grad))
        print("")
        optimizer.update()

    return model.l1.W.data[0], model.l1.b.data[0]

nxs, nys, xmeans, xstds, ymean, ystd = normalize(xs, ys)
nW, nb = train(nxs, nys)
W, b = denormalize_params(nW, nb, xmeans, xstds, ymean, ystd)

print("result of training W = %s, b = %f" % (np.array_str(W), b))
theta = normal_equation(xs, ys)
print("result of normal equation W = [%f, %f], b = %f" % (theta[1][0], theta[2][0], theta[0][0]))

実行結果

xmeans = [ 1.75900006  0.78199995]
xstds = [ 0.07020684  0.10979983]
ymean = 0.735000
ystd = 0.134257
=== Epoch 1 ===
loss = 0.644343
model.l1.W.data = [[ 0.23022108  0.00259643]]
model.l1.W.grad = [[-1.29441714 -1.51606596]]
model.l1.b.data = [ 0.]
model.l1.b.grad = [ -8.94069672e-08]

=== Epoch 2 ===
loss = 0.234075
model.l1.W.data = [[ 0.87742966  0.76062942]]
model.l1.W.grad = [[ 0.90113986  0.76939297]]
model.l1.b.data = [  4.47034836e-08]
model.l1.b.grad = [ -1.19209290e-07]

...

=== Epoch 100 ===
loss = 0.009937
model.l1.W.data = [[ 0.53428751  0.57989436]]
model.l1.W.grad = [[  6.89178705e-08   1.86264515e-08]]
model.l1.b.data = [ -8.42846930e-08]
model.l1.b.grad = [ -5.40167093e-08]

result of training W = [ 1.02172303  0.70906305], b = -1.616698
result of normal equation W = [1.021702, 0.709063], b = -1.616678

解説

予測式は以下のとおりになります(係数は小数点4位で四捨五入）。

y = 1.0217x_1 + 0.7091x_2 - 1.6167

result of training が最急降下法による推定、result of normal equation が正規方程式による推定となっています。

感想

徒然なる感想を少し。

重回帰式を chainer.links.Linear クラスで表現しています。これは Link クラスから派生したクラスです。

上のコードではごちゃごちゃと正規化(Normalization)をやっているのがわかるかと思います。本当は、複雑になるのでやりたくなかったのですが、これをやらないと学習が遅々として進まない（10万回もループを回す必要がある）ためです。正規化(Normalization)の重要さを嫌というほど思い知らされました。

私は、Python もあまりわかっていない状態なのですが、今回のコードを書くのに、Chainer のソースコードを読み込みました。Variable クラスと Function クラスにおける backward メソッドの挙動の概略を理解できたのは、大きな収穫でした。Chainer のソースコードはとても簡潔で読んでいて楽しいです。

Python では、インデントでコードのブロックが一目瞭然なのはよいですね。この点は、Haml に出会ってから大きく考えが変わりました。Ruby を愛する私ですが、ひょっとしたら（やはり、というべきか）、Ruby より Python のコードの方が読みやすいかもしれない…とか思ってしまいました。

参考文献

Chainer 本家ドキュメント