改訂履歴
2021.09.13 タイトル変更
#１．今回の目的
はじめまして。Chimoと申します。
今回Qiitaに初めて記事を書いてみようと思います。
会社での調べものやE資格の勉強の際はいつもQiitaを除いており、皆様の記事に大変お世話になりました。
私ごときが投稿するのは非常に緊張しますが、頑張ります。
記事に関してご意見やご感想ございましたら、勉強になりますのでコメントいただけますと幸いです。
さて、今回は記事を書く練習を兼ね、Boston Housingのデータセットを用いて
住宅価格の線形回帰モデルを作ることを目標にしたいと思います。
#２．データセット読み込み「Boston Housing」
scikit-learnに付属するデータセットから、Boston Housingを読み込みます。
import numpy as np
import matplotlib.pyplot as plt
#%matplotlib inline で アウトプット行にグラフを表示できる様になるらしい
%matplotlib inline
from sklearn.datasets import load_boston
以下の様に変数XとYに読み込ませます。
boston = load_boston()
X=boston.data #各特徴量
Y=boston.target #目的変数
fname=boston.feature_names #各特徴量の名称
X,Y,fnameのデータを取り出すと、ちゃんとデータが入っていることが確認できます。
X[:,1] #スライスで0行目を表示
(結果)
array([ 18. , 0. , 0. , 0. , 0. , 0. , 12.5, 12.5, 12.5,
12.5, 12.5, 12.5, 12.5, 0. , 0. , 0. , 0. , 0. ,
~中略~
0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
0. , 0. ])
Y[:5] #0~4個目までの要素を表示
(結果)
array([24. , 21.6, 34.7, 33.4, 36.2])
fname #各特徴量の名前
(結果)
array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
'TAX', 'PTRATIO', 'B', 'LSTAT'], dtype='<U7')
#３．関数作成
##３．１ 散布図表示
以下の通り散布図を表示させるための関数を作成します。
def plotdata(X,Y,num):
plt.xlabel(fname[num],fontsize=14)
plt.ylabel('住宅価格[千ドル]',fontsize=14,fontname="MS Gothic")
plt.title(fname[num] + 'と住宅価格[千ドル]の関係',fontsize=14,fontname="MS Gothic")
plt.scatter(X,Y,c='b')
plt.draw()
return
作った関数で"ZN"(*0)と価格の散布図を表示させると以下の通りです。
plotdata(X[:,1],Y,1)
(結果)
(*0)ZNは25,000 平方フィート以上の住居区画の占める割合
##３．２ 仮説
次の式で仮説:h(線形回帰モデル)を定義します。
h = X * a + b
ここでXは入力データ，a,bはパラメータを示します。
def h(X,a,b):
h = X * a + b
return h
##３．３ 目的関数
以下の通り目的関数:Jを実装します。(数式の載せ方がわからない・・ｗ)
def J(h,Y):
m = Y.size
J = np.sum(( h - y )**2) / (2*m)
return J
#４．線形回帰の学習
##４．１ 用いる特徴量
13個の特徴量を散布図で見ると"RM"(*1)が住宅価格と相関してそうなので，
適当ですが特徴量はこちらを選ぶことにします。
plotdata(X[:,5],Y,5)
(結果)
(*1)RMは住居の平均部屋数
##４．２ 学習の実施
３で作成した仮説と目的関数の関数を用いて，最急降下法によりパラメータa,bを学習させます。
ここで再急降下法は目的関数Jを各パラメータで偏微分して求めます(数式の載せ方がry)。
結果から，学習したパラメータa,bと，コスト関数の下がり具合を示したグラフが得られます。
グラフによると，100~125回くらい学習すれば十分な様ですね。
m = X.size
c = [] #costの累計に用いる配列をクリア
alpha = 0.01 #学習率
a = 0 #パラメータを初期化
b = 0
for i in range(200): #200回学習
H = h(X[:,5], a, b) #仮説にRM，パラメータ初期値を代入
J = j(H,Y) #目的関数の値を計算
c.append( J )#目的関数の値を累計
a = a - alpha * np.sum(( H - Y) * X[:,5]) / m # 最急降下法によりパラメータを学習
b = b - alpha * np.sum( H - Y) / m
plt.plot(c[0:200])#目的関数の値の類型をプロット
print('パラメータa:',a)
print('パラメータb:',b)
(結果)
パラメータa: 3.567989796463954
パラメータb: 0.4941873961329116
#５．学習結果のプロットと考察
##５．１ 学習結果のプロット
４．２にて学習した仮説h_resultをRMの散布図上に重ねてみます。
h_result = a * X[:,5] + b
plotdata(X[:,5],Y,5)
plt.xlim([0,10])
plt.ylim([0,52])
plt.plot(X[:,5],h_result, c='red')
(結果)
学習回数を10→20→・・・90→100→200(おまけ）に連れ、仮説h_resultは以下gifの様にうごきました。
##５．２ 考察
RMのプロットに対し、感覚的にはh_resultの傾きはもう少しきつい方がよさそうに思えます。
学習データの前処理を行わずに学習させたので、外れ値の影響を受けているものと推測しています。
#６．感想
初めて記事を書いてみましたが、とても楽しかったです。
数式の載せ方については調べておきたく思います・・・
次回以降でRMの前処理後の回帰や住宅価格の予測，重回帰等にも挑戦してみたいと思います。
以 上