scikit-learnで多項式回帰する方法

前回の記事で書いたのがNumPyで多項式回帰する方法だったので、今回はscikit-learnで行う方法を紹介します。
参考: NumPyで多項式回帰

NumPyの方法は、1変数の多項式に特化していたので、変数が1個しかない場合は非常に手軽に使えました。
ただ、実際は $x_0$ の多項式に加えて $x_1$, $x_2$ の変数も使って回帰したいとか、
$x_0, x_1$ を組み合わせた $x_0 * x_1$ みたいな項も入れたいとかいろんなケースがあると思います。
そのような場合は、scikit-learnの利用が検討できます。

と言っても、scikit-learnに多項式関数のモデルが実装されているわけではなく、
実際は前処理だけやってくれるモデルと、通常の線形回帰のモデルを組み合わせて使うことになります。
(このめんどくささが、1変数ならNumPyを推す所以です。)

多項式の特徴量生成には、 sklearn.preprocessing.PolynomialFeatures を使います。

試しに、3変数のデータ4セットに対して、2次までの項を生成してみたコードが次です。


import numpy as np
from sklearn.preprocessing import PolynomialFeatures

X = np.arange(12).reshape(4, 3)
print(X)
"""
[[ 0  1  2]
 [ 3  4  5]
 [ 6  7  8]
 [ 9 10 11]]
"""
poly = PolynomialFeatures(2)  # 次数2を指定
X_poly = poly.fit_transform(X)
print(X_poly)
"""
[[  1.   0.   1.   2.   0.   0.   0.   1.   2.   4.]
 [  1.   3.   4.   5.   9.  12.  15.  16.  20.  25.]
 [  1.   6.   7.   8.  36.  42.  48.  49.  56.  64.]
 [  1.   9.  10.  11.  81.  90.  99. 100. 110. 121.]]
 """

元のデータを$x_0, x_1, x_2$ とすると、
定数$1$,$x_0, x_1, x_2, x_0^2, x_0x_1, x_0x_2, x_1^2, x_1x_2, x_2^2$ のデータが生成されているのがわかります。

何番目のデータがどういう演算で生成された項なのか、という情報は、powers_ という属性に保有されています。


print(poly.powers_)
"""
[[0 0 0]
 [1 0 0]
 [0 1 0]
 [0 0 1]
 [2 0 0]
 [1 1 0]
 [1 0 1]
 [0 2 0]
 [0 1 1]
 [0 0 2]]
"""

定数1の項はいらないな、という時は、 include_biasにFalseを指定して、
PolynomialFeatures(2, include_bias=False)とすれば出てきません。

あとは、この生成されたデータを使って回帰分析を行えば、 scikit-learn を用いた多項式回帰の完成です。

NumPyで多項式回帰

前回の記事で、NumPyの多項式オブジェクトを紹介したので、ついでに多項式回帰を行う方法を紹介したいと思います。
1変数の多項式回帰に関してはscikit-learnを使うよりもNumPyの方が手軽なのでおすすめです。

使うのは、 numpy.polyfit というメソッドです。
これの引数に、回帰したいデータセットのx座標とy座標をそれぞれリストで渡し、3つ目の引数で回帰する次元を渡すだけです。
戻り値は回帰した結果の係数が配列で得られます。


import numpy as np

# 真の関数
f = np.poly1d([2/9, -3, 9, 0])

# ノイズを加えて10点サンプルを取得する
x_sample = np.array(range(10))
np.random.seed(1)
y_sample = f(x_sample) + np.random.randn(10)

# 3次多項式で回帰した係数を取得する
c = np.polyfit(x_sample, y_sample, 3)
print(c)
# [ 0.19640272 -2.60339584  7.33961579  1.29981883]

簡単ですね。

この戻り値ですが、高次の項の係数から順番に格納されており、前回の記事で紹介した多項式オブジェクトを作成するnumpy.poly1dの引数としてそのまま渡すことができます。
とても便利です。

これを使って、回帰した多項式も含めてグラフにプロットしてみます。
機械学習のテキストによくある、次数が低くて学習できてないパターンと、高くて過学習してるパターンを出してみました。


import matplotlib.pyplot as plt

# プロット用のxメモリ
x = np.linspace(-0, 9, 101)
y = f(x)
fig = plt.figure(figsize=(10, 10), facecolor="w")

for i, d in enumerate([0, 1, 3, 9], 1):
    # d次関数で回帰した係数
    c = np.polyfit(x_sample, y_sample, d)
    # d次関数オブジェクトに変換
    g = np.poly1d(c)
    ax = fig.add_subplot(2, 2, i)
    ax.plot(x, y, label="真の関数")
    ax.plot(x, g(x), label=f"d={d}")
    ax.scatter(x_sample, y_sample, label="標本")
    ax.legend()

plt.show()

出力がこちら。

どこかでみたことあるような図がバッチリ出ました。
真の関数の次数である3が一番もっともらしくフィットしていますね。

NumPyの1変数多項式クラス

Pythonで多項式を扱う時はSympyを使うものだと思い込んでいたのですが、
NumPyにも多項式クラスが用意されいるのを見つけ、しかもかなり使い勝手が良かったので紹介します。

ドキュメント: numpy.poly1d

これを使えば簡単に多項式を生成し、多項式間の演算を行ったり値を代入したりすることができます。
まず、多項式の生成は次の2種類のやりかでできます。

– 高次の項から順番に係数を指定する。
– その多項式の根を指定する。(この場合、再高次の項の係数は1)。

$2x^3-4x^2-22x+24$ と $(z-1)(z+2)=z^2+z-2$ をそれぞれの方法で定義したのが次のコードです。


import numpy as np

# 係数を指定する場合は、高次の項から順番に指定する
p1 = np.poly1d([2, -4, -22, 24])
print(p1)
"""
   3     2
2 x - 4 x - 22 x + 24
"""

# r=True を指定することで、根を指定して生成できる
# variable で変数の文字も指定できる(デフォルトはx)
p2 = np.poly1d([-2, 1], r=True, variable="z")
print(p2)
"""
   2
1 z + 1 z - 2
"""

係数のリストはc、次数はo、根の一覧はr、という属性でそれぞれアクセスすることができます。
係数と根はその多項式オブジェクトを作る時に指定した方法に関係なく取得できて便利です。
方程式を解くのもこれでできますね。


# 係数のリスト
print(p1.c)
# [  2  -4 -22  24]

# 次数
print(p1.o)
# 3

# 根
print(p1.r)
# [-3.  4.  1.]

また、指定した次数の係数は辞書と同じように[次数]でアクセスできます。


# x^2 の係数
print(p1[2])
# -4

値の代入は通常の関数と同じように、 多項式オブジェクト(代入したい値) で計算できます。


# p1 に 2 を代入
print(p1(2))
# -20

このほか、 通常の +, -, * の 演算子で演算もできます。
/ は割り算ですが、商と余りを返してくれます。とても便利ですね。


p = np.poly1d([2, 7, 1, 4, 3])
q = np.poly1d([1, 3, 5])

print(p+q)
"""
   4     3     2
2 x + 7 x + 2 x + 7 x + 8
"""

print(p-q)
"""
   4     3     2
2 x + 7 x + 2 x + 7 x + 8
"""

print(p*q)
"""
   6      5      4      3      2
2 x + 13 x + 32 x + 42 x + 20 x + 29 x + 15
"""

print(p/q)
# (poly1d([  2.,   1., -12.]), poly1d([35., 63.]))

このほかさらに、polyderとpolyintで微分と不定積分も用意されています。


# 微分
print(np.polyder(p))
"""
   3      2
8 x + 21 x + 2 x + 4
"""

# 不定積分
print(np.polyint(p))
"""
     5        4          3     2
0.4 x + 1.75 x + 0.3333 x + 2 x + 3 x
"""

例は省略しますが、2個目の引数に自然数を渡せばn回微分や、n回積分もやってくれます。
積分の方は、3つ目の引数に積分定数を渡すこともできます。

数値をゼロ埋めして桁数を揃える

桁数の揃ったIDを振るときや、表の見栄えを整える時など、桁数が少ない数字の左側に0をくっつけて表示することがあります。
そうそう頻繁にあることではないので、これまでそういう操作が必要な時は、
単純に0をたくさんくっつけて規定文字数を右側から取り出す関数を作って対応していました。

例えば次のようなメソッドを定義していました。
例として 123 を 0埋めして6桁にしています。


def zero_padding(n, length):
    m = "0" * length + str(n)
    return m[-length:]


print(zero_padding(123, 6))
# 000123

普段の用途だとこれでもあまり困らないのですが、実はPythonにはゼロ埋め専用の関数が用意されていたのに気づいたのでそちらを紹介します。
文字列オブジェクトに定義されている、zfillメソッドがそれです。

ドキュメント: str.zfill(width)
str型が持ってるメソッドなので、数値型のデータに適用するにはstrに変換してから呼び出す必要があります。

このメソッドは上のコードで僕が定義した単純な関数に比べて次の2点で優れています。
– 数字がwidth桁未満の時は、元の数字をそのまま返す。(上の方の桁をロストしない)
– 文字列の先頭が+か-の符号の場合、符号を先頭としてその後ろを0埋めする。
対応する符号は+と-だけで、±はダメみたいでした。

挙動を確認するため、いくつかのパターンでやってみます。


# 123をゼロ埋めして7桁にする
print("123".zfill(7))
# 0000123

# 幅が元の数値の桁数未満なら元の数値をそのまま返す。
print("12345".zfill(3))
# 12345

# 先頭が符号の場合は符号と数字の間をゼロ埋めする。符号も結果の文字数にカウントされる。
print("-5678".zfill(8))
# -0005678

# 小数点も文字数に数えられる。
print("+12.34".zfill(6))
# +12.34

# 数字以外の文字列に対しても使える。
print("abc".zfill(5))
# 00abc

小数や負の数をゼロパディングする機会にはまだ遭遇したことないのですが、覚えておいて損はないと思います。

Matplotlibでヒストグラムを描く時に各binのレンジを明示的に指定する方法

そもそも、Tableauなどを使えばこんな手間もないのですが、Python(Matplotlib)でヒストグラムを描く時に、各ビンの区間を指定したいことがよくあります。
0かはじめて0.5区切りにしたいとか20区切りにしたいとかです。
Matplotlibのhistメソッドでは、bins引数で、binの引数を指定でき、range引数でヒストグラムに描写する幅を指定できるので、
僕はこれまではこの二つを組み合わせて使うことで、想定通りのヒストグラムを描いていました。

試しに、 0 〜 300のデータを 20区切りで、15本のbinでヒストグラムに表示するコードがこれです。
hist メソッドの戻り値で binの区切り位置が取れるので、そちらを確認し、出力の図は省略します。
参考: matplotlibのhist()の戻り値


import matplotlib.pyplot as plt
from scipy.stats import beta

# データ生成
beta_frozen = beta(a=1, b=1, scale=300)
data = beta_frozen.rvs(100)

fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
ns, bins, _ = ax.hist(data, bins=15, range=(0, 300))
print(bins)
# [  0.  20.  40.  60.  80. 100. 120. 140. 160. 180. 200. 220. 240. 260. 280. 300.]

これ、もしreangeを指定せずに、bins=15だけ指定しているととても中途半端なところで区切られます。


fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
ns, bins, _ = ax.hist(data, bins=15)
print(bins)
"""
[  0.59033775  20.53472376  40.47910977  60.42349579  80.3678818
 100.31226782 120.25665383 140.20103985 160.14542586 180.08981188
 200.03419789 219.9785839  239.92296992 259.86735593 279.81174195
 299.75612796]
"""

さて、上記のようなblog記事ように自分で生成したデータなど取りうるレンジがわかりきってるものであれば、上記のやり方でも問題なのですが、実データでは少しだけ面倒です。
レンジとデータ量を確認して、何本くらいのbinを指定すれば切りの良い区切りで可視化でできるか考えないといけません。

しかしドキュメントあたらめて読んでみると、bins引数で、ビンの本数ではなく、区間を配列で指定できることがわかりました。
参考: matplotlib.pyplot.hist

bins : int or sequence or str, optional

If an integer is given, bins + 1 bin edges are calculated and returned, consistent with numpy.histogram.
If bins is a sequence, gives bin edges, including left edge of first bin and right edge of last bin. In this case, bins is returned unmodified.

これを使うと、例えば20区切りで可視化したい、と言った時は次のような書き方ができます。


# 300を含めるため、2個目の引数は301にしました。
bins = range(0, 301, 20)
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
ns, bins, _ = ax.hist(data, bins=bins)
print(bins)
# [  0  20  40  60  80 100 120 140 160 180 200 220 240 260 280 300]

楽ですね。
注意点として、 bisで渡した配列の区間の外側のデータは可視化されないということがあります。
bins = range(0, 300, 20) とすると、binsの配列は、
[ 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280]
になるので、 280~300のデータは可視化されません。
区間の下側も同様です。
なので、binsを指定する時に、データが全部その範囲に含まれているのかは確認しておく必要があります。

ちなみに豆知識ですが、各区間は最後(最大)の区間以外は、左の値を含み右の値を含みません、
[20, 40) の区間であれば、 $20 \leq x < 40$ のデータを数えます。
ただし、最後(最大)の区間に限って、右の値も含み、最大のbinが
[280, 300] の区間であれば、 $280 \leq x \leq 300$ のデータを数えます。

使いどこがすぐには思いつかないのですが、 bins を配列で指定する場合は、等間隔以外の区切りもサポートされているということも覚えておきましょう。

Scipyで対数正規分布

対数正規分布の話の続きです。
参考: 前回の記事

Pythonで対数正規分布を扱う場合、(もちろんスクラッチで書いてもいいのですが普通は)Scipyに実装されている、
scipy.stats.lognormを使います。
これに少し癖があり、最初は少してこずりました。

まず、対数正規分布の確率密度関数はパラメーターを二つ持ちます。
次の式の$\mu$と$\sigma$です。
$$
f(x) =\frac{1}{\sqrt{2\pi}\sigma x}\exp\left(-\frac{(\log{x}-\mu)^2}{2\sigma^2}\right).
$$

それに対して、scipy.stats.lognorm の確率密度関数(pdf)は、
s, loc, scale という3つのパラメーターを持ちます。 ややこしいですね。

正規分布(scipy.stats.norm)の場合は、 loc が $\mu$に対応し、scaleが$\sigma$に対応するので、とてもわかりやすいのですが、lognormはそうはなっていなっておらず、わかりにくくなっています。
(とはいえ、ドキュメントには正確に書いてあります。)

Scipyの対数正規分布においては、確率密度関数は
$$
f(x, s) = \frac{1}{sx\sqrt{2\pi}}\exp\left(-\frac{\log^2{x}}{2s^2}\right)
$$
と定義されています。
そして、$y=(x-loc)/scale$と変換したとき、lognorm.pdf(x, s, loc, scale) は、 lognorm.pdf(y, s) / scale と等しくなります。(とドキュメントに書いてあります。)
わかりにくいですね。

$\mu$と$\sigma$とはどのように対応しているのか気になるのですが、上の2式を見比べてわかるとおり、(そしてドキュメントにも書いてある通り、)
まず、引数のsが、パラメーターの$\sigma$に対応します。 (scaleじゃないんだ。)
そして、引数のscaleは$e^{\mu}$になります。(正規分布の流れで、locと$\mu$が関係してると予想していたのでこれも意外です。)
逆にいうと、$\mu=\log{scale}$です。

locはまだ登場していませんが、一旦ここまでの内容をプログラムで確認しておきます。
lognorm(s=0.5, scale=100) とすると、 $\sigma=0.5$で、$\mu=\log{100}\fallingdotseq 4.605$の対数正規分布になります。
乱数をたくさん取って、その対数が、期待値$4.605…$、標準偏差$0.5$の正規分布に従っていそうかどうかみてみます。


from scipy.stats import lognorm
import matplotlib.pyplot as plt
import numpy as np

data = lognorm(s=0.5, scale=100).rvs(size=10000)
log_data = np.log(data)
print("対数の平均:", log_data.mean())
# 対数の平均: 4.607122120944554
print("対数の標準偏差:", log_data.std())
# 対数の標準偏差: 0.49570170767616645

fig = plt.figure(figsize=(12, 6), facecolor="w")
ax = fig.add_subplot(1, 2, 1, title="元のデータ")
ax.hist(data, bins=100, density=True)
ax = fig.add_subplot(1, 2, 2, title="対数")
ax.hist(log_data, bins=100, density=True)
# 期待値のところに縦線引く
ax.vlines(np.log(100), 0, 0.8)

plt.show()

出力された図がこちらです。

想定通り動いてくれていますね。

残る loc ですが、 これは分布の値を左右に平行移動させるものになります。
先出の$y=(x-loc)/scale$ を $x = y*scale + loc$ と変形するとわかりやすいかもしれません、

わかりやすくするために、locに極端な値($\pm 1000$と$0$)を設定して、乱数をとってみました。
(sとscaleは先ほどと同じです。)


locs = [-1000, 0, 1000]

fig = plt.figure(figsize=(18, 6), facecolor="w")
for i, loc in enumerate(locs, start=1):
    data = lognorm(s=0.5, scale=100, loc=loc).rvs(size=1000)
    ax = fig.add_subplot(1, 3, i, title=f"loc={loc}")
    ax.hist(data, bins=100, density=True)

plt.show()

出力がこちらです。

分布の左端が それぞれ、locで指定した、 -1000, 0, 1000 付近になっているのがみて取れると思います。

対数正規分布

最近とあるデータのモデリングのために対数正規分布を使うことがありました。
そのとき使ったScipyのlognormには少し癖があったのでそれを紹介しようと思ったのですが、その前に対数正規分布そのものについて紹介させてもらいます。

対数正規分布とは一言で言ってしまうと、「それに従う確率変数の対数を取ったら正規分布に従うような確率分布」です。
正規分布の対数ではなく、対数を取ったら正規分布なので、注意が必要です。

もう少し正確に書きます。
確率変数$Y$が、正規分布$N(\mu, \sigma^2)$に従うとします。
このとき、$X=e^{Y}$は対数正規分布に従います。

定義に従って確率密度関数$f(x)$を導出しておきましょう。
まず、正規分布 $N(\mu, \sigma^2)$ の確率密度間数$g(y)$は、
$$
g(y) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(y-\mu)^2}{2\sigma^2}\right)
$$
です。
あとは$y=\log{x}$と$\frac{dy}{dx}=\frac{1}{x}$ に気をつけて、次のように計算できます。

$$
\begin{align}
f(x) &= g(y)\frac{dy}{dx}\\
&=g(\log{x})\frac{1}{x}\\
\therefore f(x) &=\frac{1}{\sqrt{2\pi}\sigma x}\exp\left(-\frac{(\log{x}-\mu)^2}{2\sigma^2}\right).
\end{align}
$$

ここで、$x$の値域は、 $0 <x < \infty$ です。
$\mu$ と $\sigma$ はパラメーターになります。
この導出の過程から明らかですが、それぞれ、対数正規分布に従う確率変数の「対数の」期待値と標準偏差です。
この二つがそのまま対数正規分布の期待値や標準偏差になるわけではないので気をつける必要があります。

対数ではない、確率変数$X$自体の期待値と分散は次の値になります。
参考: 対数正規分布 – Wikipedia
$$
\begin{align}
E[X] &= e^{\mu+\sigma^2/2}\\
V[X] &= e^{2\mu+\sigma^2}(e^{\sigma^2}-1)
\end{align}
$$

EC2の不要なインスタンスの消し方

とても単純な話なのですがかなり迷ったので備忘録的に残しておきます。

「インスタンスの削除」的なメニューがあるもんだと思ってずっと探していましたが、
実際には 「インスタンスの終了」 がそれにあたります。 (削除という文言がないのでなかなか気づきませんでした。)

消したいインスタンスを右クリックし、「インスタンスの状態」 => 「インスタンスの終了」 と選択することで消せます。
消してもしばらくはコンソール上に残るようです。

「終了保護」が設定されていると、終了ができないので、本当に消したい場合は事前に外しておきましょう。

pandasのDataFrameの欠損値をfillnaで埋める時の小技

DataFrameの欠損値(NoneやNaN)を定数で埋める時、fillnaというメソッドをよく使っていましたが、
定数で埋める以外にもいろいろ指定ができることがわかったので紹介します。

ドキュメントはこちらです。
pandas.DataFrame.fillna — pandas 1.1.2 documentation

まず、サンプルとして欠損値を含むDataFrameを作っておきます。


import pandas as pd
import numpy as np
df = pd.DataFrame(
    [[np.nan, 2, np.nan, 0],
     [3, 4, np.nan, 1],
     [np.nan, np.nan, np.nan, 5],
     [np.nan, 3, np.nan, 4]],
    columns=list('ABCD')
)
print(df)
"""
     A    B   C  D
0  NaN  2.0 NaN  0
1  3.0  4.0 NaN  1
2  NaN  NaN NaN  5
3  NaN  3.0 NaN  4
"""

fillnaの一番基本的な使い方は定数を指定するものです。
例えば、0を渡せばNaNを0に置き換えたDataFrameを返します。(inplace=True も指定すれば、データフレームそのものを書き換えます)


print(df.fillna(0))
"""
     A    B    C  D
0  0.0  2.0  0.0  0
1  3.0  4.0  0.0  1
2  0.0  0.0  0.0  5
3  0.0  3.0  0.0  4
"""

いつもはこうやって、数値なら0、文字列なら空文字列(“”)で埋めるだけの使い方をしていました。

しかし、ドキュメントを読むと、定数で埋める以外にもいろいろできます。

まず、定数ではなく、 列名: 値 の辞書を渡すことで、列ごとに別々の値で埋めることができます。


fill_values = {'A': 0, 'B': 1, 'C': 2, 'D': 3}
print(df.fillna(value=fill_values))
"""
     A    B    C  D
0  0.0  2.0  2.0  0
1  3.0  4.0  2.0  1
2  0.0  1.0  2.0  5
3  0.0  3.0  2.0  4
"""

そしてこれは、辞書の代わりに、indexに列名、valueに値を持つSeriesを渡しても同じ挙動になります。


fill_values_sr = pd.Series(fill_values)
print(fill_values_sr)
"""
A    0
B    1
C    2
D    3
dtype: int64
"""

print(df.fillna(value=fill_values_sr))
"""
     A    B    C  D
0  0.0  2.0  2.0  0
1  3.0  4.0  2.0  1
2  0.0  1.0  2.0  5
3  0.0  3.0  2.0  4
"""

これを応用すると、各列をその列の平均値や最大値、最小値で埋めることも簡単にできます。
平均値でやってみたコードが次です。


print(df.fillna(value=df.mean()))
"""
     A    B   C  D
0  3.0  2.0 NaN  0
1  3.0  4.0 NaN  1
2  3.0  3.0 NaN  5
3  3.0  3.0 NaN  4
"""

さて、ここまでは列単位である一定の値で欠損値を補完する方法でしたが、時系列データなどを扱っている時は、
一定の値ではなく、直前や直後の値で埋めたいこともあると思います。

実はこの fillna はそのようなケースにも対応しており、
method という引数に “backfill”か”bfill”を指定すれば直後の値、
“pad”, “ffill”を指定すれば直前の値で埋めることができます。

bfillとffillだけですがそれぞれ試してみたのが次のコードです。


print(df.fillna(method="bfill"))
"""
     A    B   C  D
0  3.0  2.0 NaN  0
1  3.0  4.0 NaN  1
2  NaN  3.0 NaN  5
3  NaN  3.0 NaN  4
"""

print(df.fillna(method="ffill"))
"""
     A    B   C  D
0  NaN  2.0 NaN  0
1  3.0  4.0 NaN  1
2  3.0  4.0 NaN  5
3  3.0  3.0 NaN  4
"""

以上のように、fillnaを使うと定数での補完(0埋めなど)以外にもいろんなことができることがわかりました。