Pythonで回帰分析

「わざわざPythonで？」ってところもあるかも知れませんが、とりあえず手を動かしてみようと、フィッシャー御大のあやめ(iris)のデータで回帰分析までやってみました。 参考にしたのは、「データサイエンティスト養成読本」。

データサイエンティスト養成読本 [ビッグデータ時代のビジネスを支えるデータ分析力が身につく! ] (Software Design plus)

下準備

まずは各種ライブラリーのインポートから。 回帰分析にはStatsModelsなど色々ライブラリあるみたいだが、教科書に合わせて今回はscikit-learnを選択。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import sklearn 
from pandas import Series, DataFrame

続いて、あやめ(iris)のデータのインポート

#データがスクリプトと同じフォルダにある仮定
iris = pd.read_csv("iris.csv")

データをざっくり俯瞰

data = [iris["SepalLength"],iris["SepalWidth"],iris["PetalLength"],iris["PetalWidth"],]
plt.boxplot(data)
plt.show()

本当はX軸上に項目名を表示したかったができず。 Pandasの扱いですかね、今後の宿題に。 図は左から"SepalLength"、"SepalWidth"、"PetalLength、"PetalWidth"の並び。

PetalLengthが一番分散が大きく分析として面白みがありそう。 念のため、各統計量の確認。

#項目毎の最大値
iris.max()

SepalLength          7.9
SepalWidth           4.4
PetalLength          6.9
PetalWidth           2.5
Name           virginica　#文字数の最大が返ってくる
dtype: object

#項目毎の最小値
iris.min()

SepalLength       4.3
SepalWidth          2
PetalLength         1
PetalWidth        0.1
Name           setosa　#文字数の最小が返ってくる
dtype: object

#項目毎の平均値
iris.mean()

SepalLength    5.843333
SepalWidth     3.054000
PetalLength    3.758667
PetalWidth     1.198667
dtype: float64

#項目毎の標準偏差
iris.std()

SepalLength    0.828066
SepalWidth     0.433594
PetalLength    1.764420
PetalWidth     0.763161
dtype: float64

PetalLengthのばらけ具合がどうなってるか、ヒストグラムに落としこんで確認してみる。

#PetalLengthをヒストグラムに落とし込み
plt.hist(iris["PetalLength"])
plt.xlabel("PetalLength")
plt.ylabel("Freq")
plt.show()

見事に山が二つ。

データの中身をもう少し把握したいので、品種ごとのPetalLengthの箱ひげ図を作成。

#Key"Name"でそれぞれのデータを抽出
setosa = iris[iris["Name"] == "setosa"]
versicolor = iris[iris["Name"] == "versicolor"]
virginica = iris[iris["Name"] == "virginica"]

#箱ひげ図の作成
data2 = [setosa["PetalLength"],versicolor["PetalLength"],virginica["PetalLength"],]
plt.xlabel("Name")
plt.ylabel("PetalLength")
plt.boxplot(data2)
plt.show()

やっぱり、項目名の表示がうまくいかずに挫折。 左から順に、"setosa"、"versicolor"、"virginica"。 品種ごとに値の分布が全然違うことがわかる。これが全体のヒストグラムのいびつさの原因ですね。

回帰分析を行う

今回の教科書がsetosaに注目した分析を行っていましたので、 僕のほうはvirginicaに注目することに。 まずは散布図行列でvirginiacのデータを俯瞰。

#散布図行列
pd.tools.plotting.scatter_matrix(virginica)
plt.tight_layout()
plt.show()

散布図行列からSepalLengthとPetalLengthの相関関係が強そうと判断。 「がく」が長けりゃ「花弁」も長いってことっすね。 ついでに相関係数も求めておきます。

#相関係数
virginica.corr()

          SepalLength   SepalWidth  PetalLength     PetalWidth
SepalLength     1.000000   0.457228   0.864225   0.281108
SepalWidth   0.457228  1.000000   0.401045   0.537728
PetalLength     0.864225   0.401045   1.000000   0.322108
PetalWidth   0.281108  0.537728   0.322108   1.000000

回帰分析を行う

virginicaのSepalLengthとPetalLengthの二変量について単回帰分してみる。

#単回帰分析用のライブラリーインポート
from sklearn import linear_model

#散布図を描く
x = virginica[["SepalLength"]]
y = virginica[["PetalLength"]]
plt.scatter(x,y, color = "green")


#回帰直線を散布図の中に描く
LR = linear_model.LinearRegression()
LR.fit(x,y,)
px = np.arange(x.min(), x.max(), .01)[:, np.newaxis]
py = LR.predict(px)
plt.plot(px, py, color="orange", linewidth = 2)

#図を見やすく
plt.xlabel("SepalLength")
plt.ylabel("PetalLength")

plt.show()

#統計量を求める
print(LR.coef_) #回帰係数
print(LR.intercept_) #切片
print(LR.score(x, y)) #決定係数

[[ 0.75008076]]
[ 0.61046798]
0.746884389018

きっと工夫次第でもう少しエレガントにできるんでしょうが、とりあえず今日はここまで。