線形回帰でタイタニックを試す 1

こんにちは、junkawaです。
Chromebookで始めるcoursera machine learning日記、及びkaggle日記です。
本記事では Colaboratoryで線形回帰の最急降下法を実装し、タイタニック問題を解いてみます。

courseraのmachine learningのweek2が終わりました。

学んだこと

- 線形回帰 (linear regression)という予測モデル・目的関数 (cost function / objective function)
- 線形回帰のパラメータθを求める、2通りの解法、最急降下法 (gradient discent)、正規方程式 (normal equation)
- どちらの解法を使うべきか、使う際のコツ (learning rate, feature scaling)
- 数値計算用プログラミング言語 octaveの使い方

参考

線形回帰は「連続値 (continuas data)」を予測する場合に使われる、と言われます。タイタニック問題の乗客の生死のような「分類 (classification)」を予測する場合には適用できないでしょうか?
(例えば、0以上0.5未満で死亡、0.5以上1以下で生存など)

colaboratoryで線形回帰を実装するついでに、これを試してみます。
試行錯誤の過程をそのまま記す予定です。

Colaboratoryの始め方

https://colab.research.google.com を開きます。

参考

Google Colaboratory事始め - Qiita

ColaboratoryでKaggleのデータを取得する方法

Google Colaboratory で始める Kaggler 生活(データ入手と提出編) - Qiita

データをダウンロードするために、Kaggleサイトでタイタニックコンペに参加しておきます。

また上記サイトに従って、Kaggleからkaggle.jsonをダウンロードし、Googleドライブにアップロードします。

!pip install kaggle


# download API key from google drive
## Original: https://colab.research.google.com/drive/1eufc8aNCdjHbrBhuy7M7X6BGyzAyRbrF#scrollTo=y5_288BYp6H1
## When you run for the first time, you will see a link to authenticate.

from googleapiclient.discovery import build
import io, os
from googleapiclient.http import MediaIoBaseDownload
from google.colab import auth

auth.authenticate_user()

drive_service = build('drive', 'v3')
results = drive_service.files().list(
        q="name = 'kaggle.json'", fields="files(id)").execute()
kaggle_api_key = results.get('files', [])

filename = "/content/.kaggle/kaggle.json"
os.makedirs(os.path.dirname(filename), exist_ok=True)

request = drive_service.files().get_media(fileId=kaggle_api_key[0]['id'])
fh = io.FileIO(filename, 'wb')
downloader = MediaIoBaseDownload(fh, request)
done = False
while done is False:
    status, done = downloader.next_chunk()
    print("Download %d%%." % int(status.progress() * 100))
os.chmod(filename, 600)


!kaggle competitions download -c titanic

上記をCoraboratlryに貼り付け、実行します。

これで、タイタニックコンペのデータを取得できました。

方針

- 説明変数 (independent variable) 同士が独立するように、数を減らす
  (例：客室ランク、運賃は相関が高いのでどちらか一方を使用する)
  参考: https://www.macromill.com/service/data_analysis/regression-analysis.html

- 説明変数の Feature scaling には、standardization (変数-平均)/標準偏差 を使う
  参考: https://en.wikipedia.org/wiki/Feature_scaling

- 年齢の補完には、乗船者の中央値を使う

- 1) 最急降下法を自作、2) 正規方程式を自作、3) scikit-learn の LinearRegression を使用、の順に試す。
  それぞれ別の記事にします

データの加工

import pandas as pd
train = pd.read_csv("train.csv")

ファイルの読み込みには、定番のpandasを使います。

7. 単純文 (simple statement) — Python 3.6.6 ドキュメント

pandas.read_csv — pandas 0.23.3 documentation

性別を値に変換

train.loc[train["Sex"] == "male", "Sex"] = 0
train.loc[train["Sex"] == "female", "Sex"] = 1

train.locで、"Sex"列が"male"となる行(train["Sex"] == "male")の、"Sex"列("Sex")を0としています。

pandas.DataFrame.loc — pandas 0.23.3 documentation

年齢を補完

年齢の無いデータに年齢を補完します。

train["Age"] = train["Age"].fillna(train["Age"].median())

"Age"列でN/Aとなっているセル(train["Age"].fillna)に"Age"列の中央値(train["Age"].median())をセットします。

pandas.DataFrame.median — pandas 0.23.3 documentation

pandas.DataFrame.fillna — pandas 0.23.3 documentation

使用しないデータを削除

Survived ：生存(1)、死亡(0)。目的変数 (dependent variable)
Pclass：客室ランク。今回は不使用。
Name：名前。今回は不使用。
Sex：性別。
Age：年齢。
SibSp：一緒に乗船した兄弟、配偶者の数。
Parch：一緒に乗船した親子の数。
Ticket：チケット番号。今回は不使用。
Fare：乗船料。
Cabin：部屋番号。今回は不使用。
Embarke：乗船した港。今回は不使用。

X = train.loc[:,["Sex","Age","SibSp","Parch","Fare"]]

"Sex","Age","SibSp","Parch","Fare" 列を、説明変数Xとします。
Xはm x n行列なので大文字とします。

y = train.loc[:,"Survived"]

"Survived"列を、目的変数yとします。
yは列ベクトル(m x 1行列)なので小文字とします。

説明変数の feature scaling

mu = X.mean(0)
sigma = X.std(0)
X_norm = (X - mu) / sigma

平均値mu、標準偏差sigma を求め、feature scaling後の説明変数 X_normを求めます。
mu、sigmaは 1 x 5 の行列です。X、X_normは trainのデータ数 x 5 の行列です。

X - mu で、X の各行とmuの要素同士で減算します。
また、(X - mu) / sigma で、X-muした後のXの各行とsigmaの要素同士で除算します。

pandas.DataFrame.mean — pandas 0.23.3 documentation

pandas.DataFrame.std — pandas 0.23.3 documentation

仮説関数 hypothesis

def hypothesis(X, theta):
 # h = theta0*x0 + theta1*x1 + theta2*x2 + ... thetan+xn
 # m is the number of training data.
 # n is the number of parameter.
 #
 # X = [x_0_1,x_1_1,...,x_n_1;
 #      x_0_2,x_1_2,...,x_n_2;
 #      ...                  ;
 #      x_0_m,x_1_m,...,x_n_m]
 # m x (n+1) matrix
 #
 # theta = [theta_0;
 #          theta_1;
 #          theta_2;
 #          ...    ;
 #          theta_n]
 # (n+1) x 1 vector
 #
 # h = [x_0_1*theta_0+x_1_1*theta_1+...+x_n_1*theta_n;
 #      x_0_2*theta_0+x_1_2*theta_1+...+x_n_2*theta_n;
 #      ...                                          ;
 #      x_0_m*theta_0+x_1_m*theta_1+...+x_n_m*theta_n]
 # m x 1 vector
 return X.dot(theta.values) 

実際のデータXとパラメータthetaから、yを予測する関数です。
thetaが上手く求まれば、未知のデータX'に対してy'を精度良く予測することができます。

精度良くyを予測できるようにパラメータthetaを様々に変化させることを、学習と呼びます。
精度良くyを予測できたかどうかの判定は、下記の目的関数Jで行います。

目的関数 J

def costFunction(X, y, theta):
 m = len(y)
 h = hypothesis(X, theta)
 squaredError = h.add(-y, axis=0)**2
 return np.sum(squaredError) / (2*m)

実際のデータXとパラメタthetaからyを予測した値hと、実際のyの誤差(h-y)の二乗 squaredError を求めます。

この二乗誤差を2mで割って平均を出します。(mで割るのが本来の平均の出し方ですが、gradientDiscentでの計算を楽にするために2mで割っています)

この目的関数 J が小さくなるに theta を変化させます。具体的な変化の方法は下記の 最急降下法 gradientDiscent で行います。

最急降下法 gradientDiscent

def gradientDiscent(X, y, theta, aplha, num_iters):
 J_history = np.ones(num_iters)
 m = len(y)
 for iter in range(num_iters):
   h = hypothesis(X, theta)
   g = X.T.dot(h.add(-y, axis=0).values) / m
   a = -g*alpha
   theta = theta.add(a.values, axis=0)
   J_history[iter] = costFunction(X, y, theta)
 return [theta, J_history]

目的関数Jの変数をthetaとします(X、yは定数)。
目的関数Jが最小(局所最適解)となるのは、thetaを偏微分した傾きgが0となる時です。

thetaを傾きgだけ変化させていくと、いずれ傾きgが0となります。
この時、thetaの変化の幅を抑えるために、傾きgに learning rate alpha を掛けた値 a だけthetaを変化させるようにします。

alphaが大きいとthetaの変化の幅が大きくなり、目的関数Jが収束しないで、thetaの変化に応じてどんどん値が大きくなり発散する可能性があります。

逆に、aplhaが小さすぎるとthetaの変化の幅が小さくなり、なかなか目的関数Jの値が変化しないで収束しない可能性があります。この場合、収束するためにnum_itersを大きくする(計算時間がかかる)必要があります。

適切なalphaを選択すると、すぐに目的関数Jが収束し、素早く目的のthetaを得ることができます。

目的関数Jが収束したことを確認するために、J_historyにイテレーション毎のJの値を保存しています。

目的関数Jの変化

def plotJ(J_history):
 x = np.arange(J_history.shape[0])
 y = J_history
 plt.xlabel("Number of iterations")
 plt.ylabel("Cost J")
 plt.plot(x,y)

X_norm = X_norm.assign(One=1).loc[:,["One","Sex","Age","SibSp","Parch","Fare"]]
theta = pd.DataFrame({"":[0,0,0,0,0,0]})

alpha = 0.1
num_iters = 100
[theta, J_history] = gradientDiscent(X_norm, y, theta, alpha, num_iters)

plotJ(J_history)

feature scalingしたX_normの先頭列に、全行が1となる列を追加します。
これは、theta0 * x0 の x0 に相当する列です。
theta * X の計算を行列で行うために必要となります。

X_norm.assign(One=1)で、全行が1となる列OneをX_normに追加します。
これを先頭列にするために、.loc[:,["One","Sex","Age","SibSp","Parch","Fare"]]としています。
X_normをNumPyにして計算するのが定石なのでしょうか？

alplha = 0.1、繰り返し回数num_iters = 100として、最急降下法 gradientDiscent()を実行します。
plotJ() で目的関数Jの収束を図示します。

図は、alpha = 0.1、num_iters = 100 とした場合の目的関数Jの変換です。
20回ほど繰り返した所で値が収束しているのが分かります。

テストデータに適用

def prepareData(data):
 data.loc[data["Sex"] == "male", "Sex"] = 0
 data.loc[data["Sex"] == "female", "Sex"] = 1
 data["Age"] = data["Age"].fillna(data["Age"].median())

データの加工を prepareData 関数にまとめます。

test = pd.read_csv("test.csv")
prepareData(test)

Xdash = test.loc[:,["Sex","Age","SibSp","Parch","Fare"]]
Xdash_norm = (Xdash - mu) / sigma
Xdash_norm = Xdash_norm.assign(One=1).loc[:,["One","Sex","Age","SibSp","Parch","Fare"]]

predictions = hypothesis(Xdash_norm, theta)

テストデータを読み込み(pd.read_csv)、prepareData()で整形します。

説明変数 Xdashから、トレーニングデータから求めた平均mu、標準偏差sigmaで feature scaling した Xdash_normを求めます。

この Xdash_normと、トレーニングデータで学習した theta を用いて、テストデータの生存率をhypothesis()で予測します。

predictions.columns = ["Survived"]
predictions[predictions >= 0.5] = 1
predictions[predictions < 0.5] = 0
predictions[["Survived"]] = predictions[["Survived"]].fillna(0.0).astype(int)

予測した列の名前をSurvivedとします。
線形回帰を使って生存値を求めたので、値は連続値となっています。
この連続値を生存(1)、死亡(0)に割り当てます。
ここでは、0.5以上を生存、0.5未満を0とします。
このままだと、生存が1.0、死亡が0.0となり、Kaggleの提出要件を満たさないので、.astype(int)で1、0にします。
値が存在しない(N/A)の行があるので.fillna(0.0)で0.0を埋めています。

output = test.loc[:,["PassengerId"]].join(predictions)
output.to_csv("submit.csv", index=False)

提出ファイル submit.csv を出力します(to_csv)。index列を削除するためにindex=Falseとしています。

Kaggleの提出要件が、PassengerId, Survivedなので、testファイルからPassengerIdを取り込んでいます。

Kaggleに提出

!kaggle competitions submit -c titanic -f submit.csv -m ""

スコアは0.77033でした。
前回のスコアが0.75598だったので、結構上昇しました。

まとめ

pythonの行列計算に慣れていないので、本筋とは関係ない所でつまづきました。Andrew Ng先生が初学者にoctaveをすすめる理由が分かりました。

DataFrameで扱うべき所とNumPyにすべき部分がまだ曖昧です。

DataFrameのコピーなどでメモリの無駄遣いや計算速度が遅い処理があるような気がしています。

次回は、正規方程式を実装して theta を求めてみます。