Chainer で重回帰

# 概要

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

# 重回帰

[この記事](http://d.hatena.ne.jp/mtxx/20100409/1270807079)からデータをお借りしました。ただ、単位を直しました。

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

| 身長 $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$ を予測します。



### コード

```python: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位で四捨五入）。

```math
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 本家ドキュメント](http://docs.chainer.org/)