生TensorFlow七転八倒記(2)：線形回帰を無意味に勾配法でやる - 渋谷駅前で働くデータサイエンティストのブログ

前回はロジスティック回帰をやったんですが、そう言えば普通の線形回帰やってなかったのでやっておきます。線形回帰は漫然とOLSでやるなら普通に逆行列計算しておしまいなんですが、それだと面白くないのであえて勾配法でやるという間抜けなことをやろうと思います。

TensorFlowでやってみる

冒頭のところはいつも通りです。今回使うのはBoston House Pricingデータセットを、学習データ406行とテストデータ100行に分けたものです。これまたいつも通りGitHubに転がしてあります。

まず色々必要なものをimportしておきます。あれ、csv入れっぱなの要らない気がする（Python普段使わない並みの感想）。。。データセットもpandasで読み込んでおきます。

import tensorflow as tf
import numpy as np
import csv
import pandas as pd
d_train = pd.read_csv("house_train.csv", sep=',')
d_test = pd.read_csv("house_test.csv", sep=',')

次に特徴量と目的変数とに分けるのですが、コメントアウトしたやり方だとダメだということがsess.run()して分かりました。もしかしたら.iloc()メソッドで1列だけ取ってくると行数の情報が取れないのかもしれません（TFが「行数分かんねーよ」というエラーを吐くので）。

# train_X = d_train.iloc[:, 0:13]
# train_Y = d_train.iloc[:, 13]
# test_X = d_test.iloc[:, 0:13]
# test_Y = d_test.iloc[:, 13]
train_X = d_train.iloc[:, 0:13]
train_Y = d_train[[13]]
test_X = d_test.iloc[:, 0:13]
test_Y = d_test[[13]]

あとはいつも通りplaceholderを入れて、変数スコープをreuse可に設定します。

x = tf.placeholder(tf.float32, [None, 13])
y = tf.placeholder(tf.float32, [None, 1])
W = tf.Variable(tf.zeros([13,1]))
b = tf.Variable(tf.zeros([1]))
tf.get_variable_scope().reuse_variables()

いよいよ線形回帰の本番。やることは前回同様、線形回帰子をmatmulで作るだけです。ただし今回は連続値での回帰なので、コスト関数はMSEです。前回に懲りて今回は組み込み済みのメソッドtf.losses.mean_squared_errorを使っています。

y_reg = tf.matmul(x, W) + b
cost = tf.losses.mean_squared_error(labels = y, predictions = y_reg)
optimizer = tf.train.GradientDescentOptimizer(0.000001).minimize(cost)

あとは愚直に回すだけです。。。なんですが、線形回帰を愚直に勾配法で回すと学習係数に対して鋭敏過ぎて、すぐ回帰係数が発散しますorz　微妙に色々試行錯誤した結果として、以下の設定になりました。

sess = tf.Session()
init = tf.initialize_all_variables()
sess.run(init)
num_epochs = 50000
for i in range(num_epochs):
  sess.run(optimizer, feed_dict = {x:train_X, y:train_Y})
print sess.run(W), sess.run(b)

[[-0.0972246 ]
 [ 0.09580392]
 [-0.10188288]
 [ 0.0308181 ]
 [ 0.02437116]
 [ 0.63828975]
 [ 0.11322235]
 [ 0.03513645]
 [ 0.08011973]
 [ 0.0013749 ]
 [ 0.38373682]
 [ 0.03405719]
 [-0.78699148]] [ 0.05819748]

あとはテストデータへの予測値を算出し、評価指標としてRMSEを出すだけです。

pred = sess.run(y_reg, feed_dict = {x: test_X})
rmse = ((pred - test_Y) ** 2).mean() ** .5
print rmse

MEDV    6.092334
dtype: float64

RMSE 6.1ぐらいという結果になりました。

Rで検算してみる

使い慣れたRだとこんな感じになります。

> d_train <- read.csv('house_train.csv')
> d_test <- read.csv('house_test.csv')
> d_train.lm <- lm(MEDV~., d_train)
> summary(d_train.lm)

Call:
lm(formula = MEDV ~ ., data = d_train)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.1571  -2.7694  -0.5041   1.8472  26.4347 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  35.644620   5.689105   6.265 9.80e-10 ***
CRIM         -0.106363   0.036058  -2.950  0.00337 ** 
ZN            0.029847   0.015479   1.928  0.05455 .  
INDUS         0.033390   0.069731   0.479  0.63232    
CHAS          2.143036   1.063504   2.015  0.04458 *  
NOX         -19.943015   4.450883  -4.481 9.78e-06 ***
RM            4.205103   0.463203   9.078  < 2e-16 ***
AGE           0.003266   0.014810   0.221  0.82558    
DIS          -1.378648   0.225020  -6.127 2.18e-09 ***
RAD           0.294674   0.073636   4.002 7.52e-05 ***
TAX          -0.011180   0.004142  -2.699  0.00725 ** 
PTRATIO      -0.966108   0.145893  -6.622 1.17e-10 ***
B             0.007136   0.003166   2.254  0.02474 *  
LSTAT        -0.539277   0.057079  -9.448  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.797 on 392 degrees of freedom
Multiple R-squared:  0.748,	Adjusted R-squared:  0.7396 
F-statistic: 89.51 on 13 and 392 DF,  p-value: < 2.2e-16

> rmse <- sqrt(sum(((pred-d_test$MEDV)^2)/nrow(d_test)))
> rmse
[1] 4.651266

あれー、Rの方がRMSE小さいぞorz　多分何か改善する余地があるんじゃないかと思うんですが、今回わざわざ勾配法でやったのはDNNで回帰する際にも同じ手順をたどることを念頭に置いてのことなので、余裕があったら検証することにします。。。

追記

「正規化した方が良いのでは」というコメントを直にいただきましたので、やってみました。

import tensorflow as tf
import numpy as np
import csv
import pandas as pd
d_train = pd.read_csv("house_train.csv", sep=',')
d_test = pd.read_csv("house_test.csv", sep=',')

train_X = d_train.iloc[:, 0:13]
train_Y = d_train[[13]]
test_X = d_test.iloc[:, 0:13]
test_Y = d_test[[13]]

from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
train_X.iloc[:, 0:13] = scaler.fit_transform(train_X)
test_X.iloc[:, 0:13] = scaler.fit_transform(test_X)

x = tf.placeholder(tf.float32, [None, 13])
y = tf.placeholder(tf.float32, [None, 1])
W = tf.Variable(tf.zeros([13,1]))
b = tf.Variable(tf.zeros([1]))
tf.get_variable_scope().reuse_variables()

y_reg = tf.matmul(x, W) + b
cost = tf.losses.mean_squared_error(labels = y, predictions = y_reg)
optimizer = tf.train.GradientDescentOptimizer(0.1).minimize(cost)

sess = tf.Session()
init = tf.initialize_all_variables()
sess.run(init)
num_epochs = 500
for i in range(num_epochs):
  sess.run(optimizer, feed_dict = {x:train_X, y:train_Y})
print sess.run(W), sess.run(b)

[[-0.91604364]
 [ 0.69040924]
 [ 0.23047093]
 [ 0.4955678 ]
 [-2.29030442]
 [ 2.96894979]
 [ 0.09214061]
 [-2.88571072]
 [ 2.53369761]
 [-1.87183881]
 [-2.09457803]
 [ 0.62107348]
 [-3.87335563]] [ 22.75295258]

pred = sess.run(y_reg, feed_dict = {x: test_X})
rmse = ((pred - test_Y) ** 2).mean() ** .5
print rmse

MEDV    4.858383
dtype: float64

学習係数も大きいままで大丈夫でepoch数も減らせて、RMSE 4.9ぐらいまで来ました。