lazy diary

数学とかコンピュータサイエンスとか

scikit-learn@Generalized Linear Models編

scikit-learnを使ってみます。
http://scikit-learn.org/stable/supervised_learning.html#supervised-learning

データはstatsmodelsの記事と同じものを使っています。ドキュメントを見た限りだとstatsmodelsより様々なモデルがサクッと使えるような印象があります。
散々いたるところで解説された内容なので詳細は他のブログに任せます。

今回も以下を前提とします。

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
・単回帰 / 重回帰
from sklearn import linear_model
#load data
data = pd.read_csv("./anscombe.csv")
data.drop("Unnamed: 0", axis = 1, inplace = True)
X = data["x1"].reshape((-1, 1))
y = data["y1"].reshape((-1, 1))

#モデルの作成とプロット
lm = linear_model.LinearRegression(fit_intercept = True, normalize = False)
lm.fit(X, y)
seq = np.arange(X.min() - 1.0, X.max() + 1.0, 0.1)
plt.plot(seq, lm.intercept_[0] + seq * lm.coef_[0])
plt.scatter(X, y, c = "red")
plt.grid(True)

#決定係数も求められます。
lm.score(X, y) # -> 0.6665

f:id:mt19:20180913144340p:plain

重回帰のデータはRのMASSパッケージからUScrimeを.csvで出力して使います。

#load data
from sklearn import linear_model
data = pd.read_csv("./UScrime.csv")
data.drop(["Unnamed: 0"], axis = 1, inplace = True)
X = data.drop(["y"], axis = 1)
y = data["y"].reshape(-1, 1)

#モデルの作成と決定係数の確認
mul_lm = linear_model.LinearRegression(fit_intercept = True, normalize = True)
mul_lm.fit(X, y)

N = len(y)
R2 = mul_lm.score(X, y)
R2_d = 1 - ((1 - R2) * (N - 1) / (N - X.shape[1] - 1))
print(R2_d) # ->0.7169312456540442

linear_model.score()で得られるのはあくまで決定係数なので、重回帰のように説明変数が増えた場合は自由度調整済み決定係数に変える必要があります。どうやら、sklearnにはそのメソッドは用意されていないようなので自分で少し計算する必要があります。ちなみにstatsmodelsは最初からそれを計算してくれます。

import statsmodels.api as sm
X = sm.add_constant(X)
statsmodels_ml = sm.OLS(y, X).fit()
statsmodels_ml.summary()
"""
R-squared:	0.803
Adj. R-squared:	0.717
"""

確かに上で計算した値と同じ値が得られますね。

・Ridge回帰

リッジ回帰は最小二乗項に正則化項(L1ノルム)を加えることで多重共線性による解の不安定さを解消する方法です。上の重回帰をリッジ回帰に変えてみましょう。

from sklearn import linear_model
data = pd.read_csv("./UScrime.csv")
data.drop(["Unnamed: 0"], axis = 1, inplace = True)
X = data.drop(["y"], axis = 1)
y = data["y"].reshape(-1, 1)
X = X - np.mean(X)
scale_y = y - np.mean(y)
scale_X = X.T #M行N列
M = scale_X.shape[0] #説明変数の数=14
N = scale_X.shape[1] #データ数=47

GCV_list = []
alpha_list = np.linspace(0, 10, num = 100) #正則化項の係数候補
for a in alpha_list:
    H = np.dot(scale_X.T, np.dot(np.linalg.inv(np.dot(scale_X, scale_X.T) + a * np.identity(M)), scale_X)) #N行N列
    temp = np.dot(np.identity(N) - H, scale_y)
    GCV = np.dot(temp.T, temp) / (np.power(1 - (np.trace(H) / N), 2) * N)
    GCV_list.append(GCV[0][0])
    
plt.plot(range(len(GCV_list)), GCV_list)
best_alpha_idx = np.argmin(GCV_list)
print(alpha_list[best_alpha_idx]) # -> 0

リッジ回帰を行なう前に、正則化項の係数をいくつにするかが問題になります。そこで『入門機械学習による異常検知(コロナ社)』では一般化交差確認法(GCV:generalized cross validation)と呼ばれる方法が提案されています。本にあるRのコードによるとその係数は1桁のとある数になるはずなのですが、私のpythonのコードでは係数0になってしまいました。これではリッジ回帰ではなく重回帰と変わりません。どこか間違っているようです。

「0になったけどこれ正しいの?」という私と似た質問を見つけたのですが係数の範囲を広げてみてはという回答が出ています。広げてみたところ係数が3桁で最小になるという信じがたい数字になってしまい如何にも怪しい。
Understanding ridge regression results - Cross Validated

掲載されているRのコードに使われているlm.Ridge()のGCVに関する部分を見たところ、先にcoefを求めていて本で説明されていた様子と少し違う気がします。
https://rdrr.io/cran/MASS/src/R/lm.ridge.R

長々と書いてしまったが、どうしようもないので素直にRIdgeCV()を使うことにします。

reg_CV = linear_model.RidgeCV().fit(X, y)
print(reg_CV.alpha_, reg_CV.coef_)
plt.scatter(range(len(y)),y - reg_CV.predict(X).reshape(-1, 1))
・Lasso回帰

こちらは正則化項としてL1ノルムを採用したものです。パラメータが0になることがありモデルの解釈がしやすいのが特徴です。

#Lasso regression
from sklearn import linear_model
reg_CV = linear_model.LassoLarsCV().fit(X, y)
print(reg_CV.alpha_) # ->1.7
plt.scatter(range(len(y)),y - reg_CV.predict(X).reshape(-1, 1)) #残差プロット

正則化項の係数を決めるにあたって、標本数に比べ説明変数の次元が多いときはLassoCV()、そうでなければLassoLarsCV()を使うのが良いみたいです。また、別の観点からの決め方として情報量規準(AICBIC)を使った方法がLassoLarsIC()のようです。

・Elastic Net

Ridge回帰とLasso回帰の両方を組み合わせたのがElastic Netです。

#Elastic Net
reg_CV = linear_model.ElasticNetCV().fit(X, y)
print(reg_CV.alpha_)
plt.scatter(range(len(y)),y - reg_CV.predict(X).reshape(-1, 1)) #残差プロット

結局のところどのモデルをとっても手順は全て同じです。あとはモデルの解釈の問題です。ちなみに自由度調整済み決定係数は、重回帰 > Ridge回帰 > Lasso回帰 > Elastic Net の順で、全て0.02ずつくらいの差でした。現在のデータに対してどれくらいフィットしているかなので何も正則化項の無い重回帰が大きくなるのは自然でしょう。


RIdge回帰Lasso回帰Elastic Netの詳細はこちらが分かりやすいと思います。
http://tekenuko.hatenablog.com/entry/2017/11/18/214317

・ロジスティック回帰

おなじみなので解説不要ですね。やってみます。データはstatsmodelsのときと同じように自分で作ったものです。

from sklearn import linear_model
#create data
x_1 = 5
y_1 = 7
x_2 = 7
y_2 = 5
x1 = x_1 + np.random.randn(50)
y1 = y_1 + np.random.randn(50)
x2 = x_2 + np.random.randn(50)
y2 = y_2 + np.random.randn(50)
#上を元にデータをpandas.dataframeで作成
x = pd.Series(np.hstack([x1, x2]))
y = pd.Series(np.hstack([y1, y2]))
c = pd.Series(np.repeat([1, 0], 50))
data = pd.DataFrame({"x":x, "y":y, "c":c})

#logistic regression
reg = linear_model.LogisticRegression(penalty = "l2", solver = "saga").fit(data[["x", "y"]], data["c"])
print(reg.intercept_, reg.coef_[0]) #->[-0.05453222] [-1.39388269  1.36970588]

#plot decision boundary
seq = np.arange(data["x"].min() - 1, data["x"].max() + 1, 0.1)
plt.scatter(data["x"][0:49], data["y"][0:49], c = "red")
plt.scatter(data["x"][50:], data["y"][50:], c = "blue")
plt.plot(seq, (-1 / reg.coef_[0][1]) * (reg.intercept_ + reg.coef_[0][0] * seq), c = "black")
plt.grid(True)
plt.xlabel("x")
plt.ylabel("y")
plt.title("logistic regression")

#accuracy
print(reg.score(data[["x", "y"]], data["c"])) #accuracyを計算するメソッドがある

f:id:mt19:20180916162821p:plain
statsmodelsのロジスティック回帰との違いとして、sklearnでは自動的に正則化項が組み込まれています。では正則化項が要らないという場合はどうしたら良いか調べてみたところ、正則化項の係数を小さくして影響を限りなく減らす方法が提案されていました。sklearnの場合デフォルト引数だとC = 1.0となっています。また係数が1 / Cとしてついていることに注意すると、Cをひたすら大きな数にすれば良いのです。
python - sklearn LogisticRegression without regularization - Stack Overflow

from sklearn import linear_model
reg = linear_model.LogisticRegression(penalty = "l2", solver = "saga", C = 100.0, max_iter = 1000).fit(data[["x", "y"]], data["c"])
reg.predict(data[["x", "y"]])

次にマルチクラスの場合を扱います。statsmodelsと違って難しいこと考える必要なしです。

#上に続いてもう1つデータ追加
x_3 = 1
y_3 = 10
x3 = x_3 + np.random.randn(50)
y3 = y_3 + np.random.randn(50)
x = pd.Series(np.hstack([x1, x2, x3]))
y = pd.Series(np.hstack([y1, y2, y3]))
c = pd.Series(np.repeat([2, 1, 0], 50))
data = pd.DataFrame({"x":x, "y":y, "c":c})

#lgistic regression
from sklearn import linear_model
reg = linear_model.LogisticRegression(solver = "newton-cg").fit(data[["x", "y"]], data["c"])

#plot decision boundary
seq = np.arange(data["x"].min() - 1, data["x"].max() + 1, 0.1)
plt.scatter(data["x"][0:49], data["y"][0:49], c = "red")
plt.scatter(data["x"][50:99], data["y"][50:99], c = "blue")
plt.scatter(data["x"][100:], data["y"][100:], c = "green")
plt.plot(seq, (-1 / reg.coef_[0][1]) * (reg.intercept_[0] + reg.coef_[0][0] * seq), c = "black")
plt.plot(seq, (-1 / reg.coef_[1][1]) * (reg.intercept_[1] + reg.coef_[1][0] * seq), c = "pink")
plt.grid(True)
plt.xlabel("x")
plt.ylabel("y")
plt.title("logistic regression")

#accuracy
print(reg.score(data[["x", "y"]], data["c"])) #->0.94

f:id:mt19:20180916165220p:plain

確率的勾配降下法

ここのテーマにふさわしいかどうかは微妙ですが、たまたまscikit-learnのGLMのテーマの中にあったので扱います。このクラスでは回帰や分類における10種類ほどの損失関数に対して確率的勾配降下法による最適化を提案してくれるようです。正直どの損失関数がどういう意味なのか知らないので私には違いは分からないのですが、とりあえず上のロジスティック回帰と同じことを別の最適化法で行なってみます。

from sklearn import linear_model
reg = linear_model.SGDClassifier(loss = "modified_huber", penalty = "l2").fit(data[["x", "y"]], data["c"])
print(reg.intercept_, reg.coef_)

#plot
seq = np.arange(data["x"].min() - 1, data["x"].max() + 1, 0.1)
plt.scatter(data["x"][0:49], data["y"][0:49], c = "red")
plt.scatter(data["x"][50:99], data["y"][50:99], c = "blue")
plt.scatter(data["x"][100:], data["y"][100:], c = "green")
plt.plot(seq, (-1 / reg.coef_[0][1]) * (reg.intercept_[0] + reg.coef_[0][0] * seq), c = "black")
plt.plot(seq, (-1 / reg.coef_[1][1]) * (reg.intercept_[1] + reg.coef_[1][0] * seq), c = "pink")
plt.grid(True)
plt.xlabel("x")
plt.ylabel("y")
plt.title("logistic regression by modified_huber")

f:id:mt19:20180917144532p:plain
問題無くできます。損失関数については要勉強ですね。

パーセプトロン

パーセプトロンです。いわゆる〇〇Net系で最もシンプルな形と見ることができるやつです。これも結局のところこれまでと方法は同じです。

from sklearn import linear_model
per = linear_model.Perceptron().fit(data[["x", "y"]], data["c"])
print(per.intercept_, per.coef_)

グラフはほとんど変わらないので略。


ここまで
http://scikit-learn.org/stable/modules/linear_model.html
にある様々なクラスを使いました。改めて見ると分かりますが、結局のところ全て手順は同じです。従って、これをデータに対して使い分けることがやはり重要かなと。そしてそのためにはそれぞれの特性を理解しておく必要があると感じました。

ここから他のクラスもまだまだ使ってみようと思っているのですが、方法が同じようなものはざっくりカットしていこうと思います。

statsmodelsと一般化線形モデル

統計の本で学んだものをstatsmodelsで実装してみました。statsmodelsの勉強も兼ねています。scikit-learnより様々な統計量の情報が得やすく、またモデルの設定も複雑にできるのが特徴かなと思います。そのぶん少し手間はかかります。
https://www.statsmodels.org/stable/index.html

以下

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import statsmodels.api as sm

を前提とします。

・単回帰 / 重回帰

データはRでdata(anscombe)として読み込んだデータを.csvで出力して使用。
このデータの詳細はwikipediaを見て下さい。
アンスコムの例 - Wikipedia

#load data
data = pd.read_csv("./anscombe.csv")
data.drop("Unnamed: 0", axis = 1, inplace = True) #行番号が入っていたので削除

for i in range(1,5):
    x_name = "x" + str(i)
    y_name = "y" + str(i)
    X = data[x_name]
    y = data[y_name]
    X = sm.add_constant(X) #定数項の追加
    model = sm.OLS(y, X) #OLSにデータを渡す
    result = model.fit() #モデルの生成
    print(result.summary(), "model" + str(i)) #modelの詳細を出力
    
    #プロット
    seq = np.arange(X[x_name].min() - 1, X[x_name].max() + 1, 0.1)
    plt.scatter(X[x_name], y)
    plt.plot(seq, seq * result.params[x_name] + result.params["const"], "r-") #r = red, - = 実線
    plt.show()

重回帰の場合はfit()でのXが行列になります。

多項式回帰

data(anscombe)の2つ目のデータは明らかに曲線上の点なので、2次で多項式回帰をしてみます。やることは単回帰と同じです。

X = data["x2"]
y = data["y2"]
#2次の項と定数項を追加します
X = pd.DataFrame(X)
X["x2_2"] = X * X
X = sm.add_constant(X)

result = sm.OLS(y, X).fit() #モデルの作成
#プロットする
seq = np.arange(X["x2"].min() - 1, X["x2"].max() + 1, 0.1)
plt.scatter(X["x2"], y)
plt.plot(seq, result.params["const"] + result.params["x2"] * seq + result.params["x2_2"] * (seq * seq), "r-")

f:id:mt19:20180822122454p:plain
かなりフィットしてしまいました(偶然です)

カーネル密度推定法(Kernel Density Estimation)

『入門機械学習による異常検知』でこの方法を知りました。statsmodelsで実装できるみたいなのでやってみます。カーネル密度推定の詳細は他のブログでお願いします。データはRからlibrary(car), Davisで読み込み.csvで出力したものを使っています。

data = pd.read_csv("./car_data.csv")
data = data[data["height"] > 100] #1つ異常値があったので削除しておきました
data = data["height"].astype(np.float) #floatにしないと後々エラーを出すのでキャストしました

#カーネル密度推定をします
kde = sm.nonparametric.KDEUnivariate(data)
kde.fit()
plt.hist(data, normed = True) #ヒストグラムを相対度数にしておきます
plt.plot(kde.support, kde.density,"r-")
plt.grid(True)
plt.show()

f:id:mt19:20180828143327p:plain
特に何も設定しない場合はカーネル関数はガウスカーネルになるらしい。他の関数を使う場合は設定は以下のようにする。

from statsmodels.nonparametric.kde import kernel_switch
list(kernel_switch.keys()) 
#使用できるカーネル関数一覧 ['biw', 'cos', 'cos2', 'tri', 'uni', 'epa', 'triw', 'gau']

kde = sm.nonparametric.KDEUnivariate(data)
kde.fit(kernel = "tri", fft = False) #ここで設定

カーネル関数を変更できました(変更後の関数がどういう関数なのかはすみませんよく知りません)
ちなみにバンド幅を自分で設定したり確認したりもできるそうです。

kde = sm.nonparametric.KDEUnivariate(data)
kde.fit(bw = 3.3) #ここで設定
kde.bw #これで確認できる...3.3
ポアソン回帰

GLMの1つとしてポアソン回帰をしました。過分散が生じていてうまくデータを表現できていませんが実装としては次のやり方でできます。
データはみどり本からdata3a.csvを利用。
http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html

import statsmodels.formula.api as smf
df = pd.read_csv("./data3a.csv")
df = pd.get_dummies(df)

#GLMの実装
df = sm.add_constant(df) #定数項の追加
X = df.drop(["y"], axis = 1)
y = df["y"]
formula = "y ~ x + f_C + f_T"
glm_model = smf.glm(formula = formula, data = df, family = sm.families.Poisson()).fit() #GLMの作成

model.summary() #これで各統計量をチェック

#プロット
plt.scatter(range(len(y)), y, c = "red")
plt.scatter(range(len(y)), model.predict(), c = "blue")

f:id:mt19:20180905105312p:plain

・ロジスティック回帰

この詳細は『一般化線形モデル入門(共立出版)』で学べます。実装は引き続きGLMを使います。背景を知らないと手続きが面倒です。二項分布なのでリンク関数はロジット関数を使うとか。もっと簡潔に書けるような気がしますが、とりあえずこれで。

#データは人工的に作ります。
x_1 = 5
y_1 = 7
x_2 = 7
y_2 = 5
x1 = x_1 + np.random.randn(50)
y1 = y_1 + np.random.randn(50)
x2 = x_2 + np.random.randn(50)
y2 = y_2 + np.random.randn(50)
plt.scatter(x1, y1, c = "r")
plt.scatter(x2, y2, c = "b")
x = pd.Series(np.hstack([x1, x2]))
y = pd.Series(np.hstack([y1, y2]))
c = pd.Series(np.repeat([1, 0], 50))
data = pd.DataFrame({"x":x, "y":y, "c":c})
data

f:id:mt19:20180905120639p:plain
このようなデータです。これの決定境界をGLMを使って引いていきます。

#ロジスティック回帰を実装する
import statsmodels.formula.api as smf
formula = "c ~ x + y"
data = sm.add_constant(data)
model = smf.glm(formula = formula, data = data, family = sm.families.Binomial()).fit()

#accuracyを計算
def f(x):
    if x >= 0.5:
        return 1
    return 0
predict = [f(x) for x in model.predict()]
accuracy = np,sum(predict == data["c"]) / len(data["c"]) #0.93でした

#決定境界を引く
#0 = intercept + x_coef * x + y_coef * yに注意
seq = np.arange(data["x"].min() - 1, data["x"].max() + 1, 0.1)
plt.plot(seq, (-1) * (model.params["x"] * seq + model.params["Intercept"]) / model.params["y"], c = "green" )
plt.scatter(x1, y1, c = "r")
plt.scatter(x2, y2, c = "b")

f:id:mt19:20180905120936p:plain
次にマルチクラスの分類をしようと思ったのですが、GLM関数を使う場合one vs allの方法以外が浮かばなかったので、statsmodelsで別の関数を探してみました。すると多項ロジットのための関数があったのでそれで実装してみます。

#3つ目のデータを作ります
x_3 = 1
y_3 = 10
x3 = x_3 + np.random.randn(50)
y3 = y_3 + np.random.randn(50)
x = pd.Series(np.hstack([x1, x2, x3]))
y = pd.Series(np.hstack([y1, y2, y3]))
plt.scatter(x1, y1, c = "r")
plt.scatter(x2, y2, c = "b")
plt.scatter(x3, y3, c = "green")

#dataframeの作成
x = pd.Series(np.hstack([x1, x2, x3]))
y = pd.Series(np.hstack([y1, y2, y3]))
c = pd.Series(np.repeat([2, 1, 0], 50))
data = pd.DataFrame({"x":x, "y":y, "c":c})

f:id:mt19:20180905151923p:plain
マルチクラス分類を行ないます。

#多項ロジットによるロジスティック回帰
import statsmodels.discrete.discrete_model as dm
multi_log = dm.MNLogit(data["c"], data[["x", "y"]]).fit()

#accuarcyを求める
predict = pd.DataFrame(multi_log.predict())

def max_index(x):
    return x.argmax()

predict_index = predict.apply(max_index, axis = 1)
accracy = np.sum(predict_index == data["c"]) / len(data["c"])
print(accuracy) #0.93でした

#決定境界を引く
x0_coeff = multi_log.params[0]["x"]
y0_coeff = multi_log.params[0]["y"]
x1_coeff = multi_log.params[1]["x"]
y1_coeff = multi_log.params[1]["y"]

#多項ロジット関数なので線形予測子とlog(p0 / p3), log(p1 / p3)であることに注意
seq = np.arange(data["x"].min() - 1, data["x"].max() + 1, 0.1)
plt.plot(seq, (-1) * (x0_coeff - x1_coeff) * seq / (y0_coeff - y1_coeff), c = "black" )
plt.plot(seq, (-1) * x1_coeff * seq / y1_coeff, c = "purple")
plt.scatter(x1, y1, c = "r")
plt.scatter(x2, y2, c = "b")
plt.scatter(x3, y3, c = "green")

f:id:mt19:20180905152159p:plain
当初は切片を用意するためにsm.add_constant()をしていたのですが、おそらく桁数の問題で計算がうまくいかなかったのでいつかまた。

・分散分析

まずは一元配置分散分析を行ないます。このための関数もstatsmodelsにはあります。データはRからdata(InsectSprays)で読み込んだものを.csvで出力して使います。

data = pd.read_csv("./insectsprays.csv")
data.drop(["Unnamed: 0"], axis = 1, inplace = True)
data["spray"].value_counts() #データの中身を調べます
"""
B    12
E    12
C    12
A    12
D    12
F    12
Name: spray, dtype: int64
"""

#可視化してみます
insect_data = data.groupby(["spray"])
colors = ["r", "g", "blue", "yellow", "black", "purple"]
cnt = 0
for values, group in insect_data:
    plt.plot(range(12), group["count"], color = colors[cnt], label = values)
    cnt += 1
plt.legend()
plt.title("insect counts")
plt.grid(True)

f:id:mt19:20180905204506p:plain
ここからが本題です。

#ANOVAに必要な線形予測子の部分を準備
import statsmodels.formula.api as smf
formula = "count ~ spray"
insect_lm = smf.ols(formula, data = data).fit()
insect_lm.summary()

#ANOVAの実装
aov_table = sm.stats.anova_lm(insect_lm, typ=2)
aov_table

これでF値が出力できます。
次に二元配置分散分析をしてみます。データはRからdata(ToothGrowth)で読み込んだものを.csvで出力しています。手順は上と同じです。

#読み込んで可視化
tooth_data = pd.read_csv("./ToothGrowth.csv")
tooth_data.drop(["Unnamed: 0"], axis = 1, inplace = True)
category = ["supp", "dose"]
groupby_data = tooth_data.groupby(category)
data_num = 10
colors = ["r", "g", "blue", "yellow", "black", "purple"]
cnt = 0

for values, group in groupby_data:
    plt.plot(range(data_num), group["len"], color = colors[cnt], label = values)
    cnt += 1
plt.legend(loc = "upper left")
plt.title("tooth length")
plt.grid(True)

f:id:mt19:20180906102318p:plain
ここからが本題です。formulaの部分を変えることで他にも様々なモデル(加法モデルなど)に適応できます。今回はフルモデルです。

#分散分析
import statsmodels.formula.api as smf
formula = "len ~ supp + dose + supp * dose"
tooth_lm = smf.ols(formula, data = tooth_data).fit()
tooth_lm.summary()

aov_table = sm.stats.anova_lm(tooth_lm)
aov_table

これでF値を得られます。

・生存時間解析(カプラン・マイヤー推定の利用)

生存時間解析と一口に言っても色々あるのですが、とりあえず『一般化線形モデル入門(共立出版)』にあるカプラン・マイヤー推定をしてみます。データはRでlibrary(MASS) data(gehan)として読み込んだものを.csvとして出力して使います。

"""
白血病患者に対する薬の効果を調べたデータ。
pair:投薬と比較対象のペア
time:生存時間
cens:0は打ち切りデータを意味する
treat:6-MPを投与したかどうか
"""
data = pd.read_csv("./gehan.csv")
data.drop(["Unnamed: 0"], axis = 1, inplace = True)

#カプラン・マイヤー推定で経験生存関数を求める
bytreat_data = data.groupby(["treat"])
num = 21
cnt = 0
color = ["red", "blue"]
for values, group in bytreat_data:
    sf = sm.SurvfuncRight(time = group["time"], status = group["cens"])
    print(values)
    print(sf.summary())
    sf.plot()
    plt.legend(loc = "upper left")
    plt.grid(True)
    plt.title("prob of the duration({})".format(values))
    cnt += 1

f:id:mt19:20180906120322p:plainf:id:mt19:20180906120336p:plain
投与したほうが元気でいられる時間が長いらしいことは分かります。切り込みは打ち切りデータであることを意味しています。本来はこのあと別のグラフからワイブル分布に近いことを察してパラメータを推定する予定だったのですが、どうやらモデリングがあまりうまくいっていなかったのでここでは省略します。

statsmodelsの特徴として述べた統計量は簡単に調べられるのですが、ここでは載せていません(summary()で調べられます)。ここで実装したことの数学的な詳細は主に『一般化線形モデル入門(共立出版)』で述べられているので興味があれば読んでみて下さい。カーネル密度推定だけは別の本です。他にもstatsmodelsには時系列データの扱いに関するメソッドが豊富に用意されているみたいのですが、あまり時系列データのことは知らないのでまた今度にします。

一般化線形モデル入門 原著第2版

一般化線形モデル入門 原著第2版

私の図書館(2018/09/18更新)

https://www.pakutaso.com/shared/img/thumb/elly20160628465420_TP_V.jpg


私が学習したものでかつ勧められそうなものを載せてあります。

Python

Python学習講座~PythonエンジニアによるPython3学習支援サイト~

Pythonの初歩から始まっています。
基礎編はいわゆる一般的な入門書にあるPythonの話で、応用編はライブラリの話です。データ分析に必要なライブラリ「matplotlib/pandas/numpy」の説明もあるので、応用編で興味のあるものをいくつか取り組むと良いと思います。
PythonエンジニアによるPython3学習サイト

・Aidemy

お手軽にPython機械学習の概観を勉強できます。
(2018/05/07追記:いつからか有料になっていました。当初は全講座が無料でした。また、当初よりコンテンツがかなり増えているようです。)
aidemy.net

・【ゼロから始めるデータ分析】 ビジネスケースで学ぶPythonデータサイエンス入門

Udemyの講座です。分析の一通りの手順を教えてくれます。
実際にSignateに出てくる問題を扱っているので、実践形式で学べるところが良し。
【ゼロから始めるデータ分析】 ビジネスケースで学ぶPythonデータサイエンス入門 | Udemy

R

・全人類がわかる統計学

Rの基礎をこのサイトで覚えました。
to-kei.net

機械学習

・Coursera:Machine Learning

この世界なら知らない人はいないであろうAndrew Ng先生の授業です。
内容は機械学習の基礎で、大半は日本語字幕があります。たまに大幅に字幕がずれていて違和感しかないものもありますが、そのときは英語で乗り切ります。WEEK1-11まであり、各WEEKの最後にはクイズやプログラミングの宿題が課されます。具体的な授業内容は別の記事に書いたので、良ければ参考に。
Coursera:Machine Learningを受講しました - lazy diary
www.coursera.org

・ゼロから作るDeep Learning――Pythonで学ぶディープラーニングの理論と実装

NNをPythonでスクラッチして理論を学んで最終的にCNNまでやってみようという本です。分かりやすいと人気の高い本ですが、第5章の計算グラフを使ったback propagationの説明は好みが分かれるような気がします。しかしながら第5章以外、特に第6章の学習に関するテクニック(最適化の方法)は知らないこともあったので勉強になりました。

統計学

統計学入門 (基礎統計学Ⅰ)

統計の赤本の名で有名な統計学の入門書です。線形代数だと行列の基本的な扱い、解析だと重積分あたりが前提知識かと思います。古い本ですがそれを感じさせない読みやすさです。分布収束といった収束の概念については深く説明しないようになっていますが、入門書なのでこれで良いと思います。日常でも聞く平均などの統計量の説明から始まり、そこに確率変数を導入して分布の説明、次に標本の話と標本からパラメータを推定したり仮説検定を行なったりする話、最後に少しですが回帰分析の話が載っています。
私はこの本を統計検定2級を受験する際の教科書として利用しました。統計検定2級の受験記も良ければ見て下さい。
統計検定2級に合格しました - lazy diary

統計学入門 (基礎統計学?)

統計学入門 (基礎統計学?)

・現代数理統計学の基礎

統計学入門』と重複する内容が多いです。しかしながら、分布収束といった収束の概念を使って細かい説明がなされています(※測度論の話は数か所で出てくる程度です)。まさに"数理"統計の本です。第8章が検定理論で一般的な統計学の本が扱うところはこのあたりまででしょう。第9章以降は応用で、「回帰モデル/リスク最適化/モンテカルロ法/ブートストラップ/確率過程」を主に扱っています。ここは説明が簡潔で難しいと感じました。他の文献を参照しながらでないと読むのが難しかったです。
最後に、この本の良い点は誤植や例題の解答、補足等を著者がWEBサイトで公開している点だと思います。数学書でこれほどありがたいことはないです。

現代数理統計学の基礎 (共立講座 数学の魅力)

現代数理統計学の基礎 (共立講座 数学の魅力)

・データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

みどり本の名で有名な統計モデリングの本です。『現代数理統計学の基礎』が検定理論を中心とした数理統計の本なら、この本はモデルというものに対する考え方を紹介した本です。「一般化線形モデル(GLM)、一般化線形混合モデル(GLMM)、マルコフ連鎖モンテカルロ法(MCMC)、階層ベイズモデル」の順番で説明がされており、なぜこのモデルではデータを表現しきれていないのか、そのために何が必要なのか、それをどう実装するのか(Rのコード付き)という流れで理解を深めることができます。本書で使われているWinBUGSは(開発が止まっているから?)古いという声もあって、下記のPyStanに移植したというありがたい記事が私は非常に参考になりました。
データ解析のための統計モデリング入門(緑本)の読書メモ(PythonとStanで)

・一般化線形モデル入門

一般化線形モデルに関する古典的名著の日本語版です。一般化線形モデルを数理統計の視点で解説しています。一般化線形モデルのパラメータ推定法を最初に説明したのち、正規線形モデル(重回帰、分散分析、共分散分析)、ロジスティック回帰(2値変数の場合からマルチクラスへ、また名義尺度、順序尺度の場合の解説もあります)、ポアソン回帰、生存時間解析、また検定理論を使ったモデル選択の話があります。最後に簡潔ではあるが時系列データ、クラスターデータにも触れています。この本を読めば一般化線形モデルの数理の側面は一通り理解できると思います。この本の内容のだいたいをstatsmodelsで実装したので良ければ参考に。
statsmodelsと一般化線形モデル - lazy diary

一般化線形モデル入門 原著第2版

一般化線形モデル入門 原著第2版


その他

・入門 機械学習による異常検知―Rによる実践ガイド

異常検知に興味が湧いたので読みました。おそらく大学で初年度に学ぶ解析学線形代数を使い慣れていないと難しいと思います。それくらいギュッと内容が凝縮された本でした。途中式はあまり行間が飛ぶことなく丁寧に書かれており、著者の親切さ(のようなもの?)をひしひしと感じました。内容は異常検知の典型的なステップをホテリング理論を追いながら学び、続いてクラスタリング、主成分分析、回帰モデル、時系列モデリングに関する異常検知モデルの解説、また途中には異常検知モデルの評価方法、そして最後には実データを扱うときにぶつかる問題についてのフォローがされています。統計学の本で学ぶ沢山の分野が詰まっています。私自身まだまだ理解しきれてない点があるので何度も読み直すこと必須のような気がしてならないです。

赤池情報量規準AICモデリング・予測・知識発見― (現在、取り組み中)

統計モデリングで頻出するAICについての本です。本書はその主役である赤池博士の業績および、博士によって切り拓かれた情報量に基づく情報処理の最先端の成果を紹介するものとして企画されたみたいです。AICを考えるようになったきっかけから始まりAICとその周辺の事柄がまとめられています。私個人としてはモデリングは本書のようにmotivationの部分の話があると理解が一気に深まるような気がしています。

赤池情報量規準AIC―モデリング・予測・知識発見

赤池情報量規準AIC―モデリング・予測・知識発見

統計検定2級に合格しました

合格したので感想を。
http://www.toukei-kentei.jp/wp-content/uploads/ToukeiAkari-300x123.png*1

統計検定とは

www.toukei-kentei.jp
そもそも統計検定って何?って方もいると思うのでザックリと説明を。

「統計検定」とは、統計に関する知識や活用力を評価する全国統一試験です。
データに基づいて客観的に判断し、科学的に問題を解決する能力は、仕事や研究をするための21世紀型スキルとして国際社会で広く認められています。
日本統計学会は、中高生・大学生・職業人を対象に、各レベルに応じて体系的に国際通用性のある統計活用能力評価システムを研究開発し、統計検定として実施します。

うむ、よく分からんということで具体的な検定内容は次のようになっています。

試験の種別 試験内容
統計検定1級 実社会の様々な分野でのデータ解析を遂行する統計専門力
統計検定準1級 統計学の活用力 ─ データサイエンスの基礎
統計検定2級 大学基礎統計学の知識と問題解決力
統計検定3級 データの分析において重要な概念を身に付け、身近な問題に活かす力
統計検定4級 データや表・グラフ、確率に関する基本的な知識と具体的な文脈の中での活用力
統計調査士 統計に関する基本的知識と利活用
専門統計調査士 調査全般に関わる高度な専門的知識と利活用手法

2/3級については、CBT方式での受験が可能です。従って、事前に申し込みをすれば指定された全国の会場でほぼ毎日受検することが出来ます。
一方で、1/準1級については受験日が決まっています。1級は11月、準1級は6月に試験が開催されるみたいです。全て合格点は60/100点だったと思います。
統計調査士、専門統計調査士についてはよく知らないので略。

この統計検定始まったのは2011年とまだ最近です。
AIやデータサイエンスの人気の高まりを受けて受験者数は増えるんでしょうかね。
かく言う私もその口なのですが。

勉強したこと

私は数学畑の人間とは言っても、お恥ずかしいことに統計学については素人です。統計??大学受験レベルのものと、大学の実験の講義で少し学んだ単回帰くらいなら分かりますという今思えばひどい知識でした。確率分布の知識はゼロです。そこで合格した人のブログをいくつか調べて、この本で学ぶことにしました。

統計学入門 (基礎統計学?)

統計学入門 (基礎統計学?)

学習まとめのブログに書評は書きましたが、とても読みやすい本でした。
この本をノートにまとめながら読むこと約2週間(※練習問題は飛ばしました)。
これだけでどれくらい点数が取れるのだろうかと過去問を解くことにしました。
日本統計学会公式認定 統計検定 2級 公式問題集[2014〜2016年]

日本統計学会公式認定 統計検定 2級 公式問題集[2014〜2016年]

ギリギリですが合格点を取ることが出来ました。実はこのとき初めて統計検定2級の範囲を見たのですが、『統計学入門』にいくつか載っていない項目があったのでそこをググりながら補うことにしました。これに1週間ほどかけました。その後、過去問の再度解き直しに1週間をかけ、安定して合格できるくらいまで知識を確認しました。

受験当日

事前に日程をこちらが決めてCBT方式で受験しました。会場に行ってパソコンの前に着席してから説明を聞いて、すぐに試験開始です。メモ用紙が貰えるでそこでガリガリ計算していきました。偶然かもしれませんが、過去問よりいくぶんか簡単でした。その場で合否が分かるので合格発表をドキドキして待つこともありません。

最後に

ここまで言及していなかったのですが、忘れてはいけないのが電卓です。この試験をサクサク乗り切るためには電卓をうまく使うことが結構重要なのではと過去問を解きながら思っていました。私は電卓といわれるとスマホの電卓しか長らく使っていませんでしたが、スマホの電卓はNGなので今回この試験のために買いました。

カシオ 電卓 時間・税計算 手帳タイプ 8桁 SL-300A-N

カシオ 電卓 時間・税計算 手帳タイプ 8桁 SL-300A-N

あまり使う機会ないですしお金が無い学生にはこれで十分でした。この試験で今まであまり使うことのなかったMC, MR, M-, M+も使わざるを得なくなりました。意外と便利ですね。


参考書を選んでしまえばあとはそれで学習して最後に過去問で演習という流れで合格できると思います。参考書は範囲をおおよそカバーしているのものの中から気に入ったのを選べば良いと思います。『統計学入門』もその1つです。多くの有名な資格試験(ex.情報処理試験)はそれに最適化された本が出版されていますが、統計検定は頼りない公式の参考書しかありません。その中で、1冊で範囲のほとんどを押さえることが出来るこの本はとてもありがたいと言えるのではないでしょうか。


次は統計検定1級に挑戦してみようと思います。頑張ります。

*1:統計検定の公式キャラクターらしいです http://www.toukei-kentei.jp/info/toukeiakari/

Coursera:Machine Learningを受講しました

機械学習を勉強している人ならどこかで聞いたことはあるであろう(たぶん?)
あの有名なAndrew Ng先生のMachine Learningを受講しました。
www.coursera.org

何を学んだのかの振り返りの意味も込めて、この記事を書いていきます。

1.受講動機

機械学習をざっくりと学びたく、当初市販の本を手にとってはみたのですがどうもしっくり来ませんでした。市販の本はたくさん数式が出てくるのは良いんですが、なぜその手法を使うのかといった所謂motivationの部分の話が無くて分かりづらかったのです。そこで本から離れて評判の良いCourseraの講義に手を出したということです。

2.講義内容

散々色んなブログでまとめられているとは思いますが、振り返りの意味も込めて自分で作ったノートを見ながら講義内容をざっくりここに箇条書きにします。

Week1

機械学習の一般的な定義と教師あり学習教師なし学習について
MATLABOctaveアルゴリズムの簡単な実装には便利
教師あり学習、回帰問題の初歩として線形回帰(単回帰)を学ぶ
最急降下法とは何か?

市販の本で分かりづらかった最急降下法がとても分かりやすく、あぁ教えるのがお上手でと思いながら動画を見てました。

Week2

・単回帰の拡張としての重回帰
・スケーリングの重要性とその実装
・単回帰/重回帰の拡張としての多項式回帰
最急降下法以外のパラメータの求め方として正規方程式の利用

なぜその手法が必要なのか、どのようなメリットがあるのかが説明されているので、すんなりと新しいことにも取り組めました。

Week3

教師あり学習、分類問題の初歩としてロジスティック回帰を学ぶ
・ロジスティック回帰にも最急降下法を適用する
・マルチクラスの分類問題への対応の仕方
過学習の恐ろしさとその対策方法(特徴量、正規化etc)
・線形回帰、ロジスティック回帰への正規化の適用

Octaveの便利さをこのあたりで痛感。

Week4

・ロジスティック回帰の欠点
ニューラルネットワークの概要とその素晴らしさ

Week5

・forward propagationとback propagation
ニューラルネットワークにおける最急降下法とエラーチェック(gradient checking)

ニューラルネットワークは線形回帰やロジスティック回帰に比べてかなり難しかったです。特にback propagationがなぜ必要なのかがなかなか理解できませんでした。Machine Learningの講義で私が一番大変だと思ったWeekですね。

Week6

・交差検定の理解
・正規化における定数λの決定方法
・学習曲線から分かること
・モデルの改善方法の決め方
・分類問題に対する機械学習のシステムデザイン(F値の重要性)

ここではアルゴリズム云々ではなく、モデルの改善にはなにが必要なのかが学べます。とても重要なんでしょうが、ニューラルネットワークほど重くないのであっさり進めました。

Week7

・サポートベクターマシーン(SVM)の実装とその良さ
非線形に対応するためのカーネルの理解
SVMによるマルチクラスの分類問題への対応
・ロジスティック回帰とSVMの使い分けについて

SVMは強力ということを何度も言われます。実際にどこかで使いたい。

Week8

教師なし学習の応用例
教師なし学習の初歩としてK-means法を学ぶ
教師なし学習におけるコスト関数とその利用
・マルチクラスにおける定数Kの決め方
・次元削減の必要性と主成分分析の概要
機械学習での主成分分析の組み込み方
・主成分の数kの選び方と組み込むときの注意

教師なし学習教師あり学習より理解は簡単でした。
次元削減は説明は簡単ですが、実際にアプリケーションに組み込むとなると大変そうです。

Week9

・異常検知の初歩とガウス分布
・異常検知アルゴリズムの評価方法と改善方法
・特徴量をガウス分布に変換する方法
・多変量ガウス分布の必要性とその実装
ガウス分布と多変量のガウス分布の使い分け
・レコメンデーションシステムの概要
協調フィルタリングの実装と実際のレコメンドの仕方
・平均標準化の利用

項目はたくさん並んでしまいましたが、ガウス分布等は統計では定番なので、統計を少し学んだ人は苦労しないと思います。レコメンデーションシステムは機械学習の定番なので知れてよかったですね。

Week10

・ここまで取り扱ってきた最急降下法(バッチ急降下法)のデメリット
・確率的急降下法とミニバッチ最急降下法の実装
・オンライン処理の特徴
・データが膨大にあるときに使われるMap Reduceについて

データをたくさん持っていることは強い(オーバーフィッティングしているとき)というのは動画内で何度も説明されます。その膨大なデータを計算コストの面からどのように扱うかが説明されます。
Andrew Ng先生が常に処理時間を気にされていることが分かります。

Week11

機械学習の応用例としてのOCRの説明
・アプリケーションのパイプラインを用意することの重要性
・現実で直面するデータ量の問題と人工データ合成の概要
機械学習エンジニアにおける時間の重要性
・シーリング分析の概要
・あなたがこの講義で学んできたことのまとめと、受講に対する感謝の言葉

機械学習を実際のアプリケーションに組み込む際の注意事項がたくさん説明されています。
最後の受講に対する感謝の言葉は涙ものですよ。

3.受講してみて

受講して本当に良かったです。とにかく分かりやすいの一言です。
学生の特権、春休みの1ヵ月使って学びましたが、勉強していて楽しかったです。動画なので本と違って行間埋める必要もありませんのでとっつきやすいです。

もし

前提条件
機械学習の基礎を学んでみたい
・英語が少し読める
歓迎条件
・行列演算、微分偏微分を知っている(これらは知らなくても動画内で簡単な説明がありますが、この先どこかで深い知識が必要になるでしょうから事前に数学書で学ぶことをオススメします)
・巷の機械学習入門書が難しすぎる

という方はオススメです。