statsmodelsでかばん検定 (自己相関の検定)

時系列データを分析をするとき、そのデータが自己相関を持つかどうかはとても重要です。
データが自己相関を持っていたらその構造を記述できるモデルを構築して、予測等に使えるからです。
逆に自己相関を持っていないと、時系列分析でできることはかなり限られます。
過去のデータが将来のデータと関係ないわけですから当然ですね。
(どちらも沖本本の 1.4 自己相関の検定から)

ということで今回行うのは自己相関の検定です。

自己相関が全てゼロという帰無仮説、つまり
$H_0:\rho_1=\rho_2=\cdots=\rho_m=0$を、
$H_1:$少なくとも1つの$k\in[1,m]$において、$\rho_k\neq0$
という対立仮説に対して検定します。

この検定はかばん検定(portmanteau test)と呼ばれているそうです。
検定量は色々考案されているそうですが、 Ljung and Box が考案されたものがメジャーとのこと。

具体的な数式や、numpyでの計算例は次回に譲るとして、とりあえずpythonのライブラリでやってみましょう。
statsmodels に acorr_ljungbox という関数が用意されています。

statsmodels.stats.diagnostic.acorr_ljungbox
(完全に余談ですが、statsmodelsのドキュメントで portmanteau という関数を探していたので、これを見つけるのに結構苦労しました)

あからさまに7点周期を持つデータを準備し、1から10までのmに対して検定を実施したコードがこちらです、
acorr_ljungbox は lb値と、p値をそれぞれのmに対して返します。


import pandas as pd
import numpy as np
from statsmodels.stats.diagnostic import acorr_ljungbox

# 7点ごとに周期性のあるデータを準備
series = pd.Series([1, 1, 1, 1, 1, 1, 5]*10)
# 乱数追加
series += np.random.randn(70)

lbvalues, pvalues = acorr_ljungbox(series, lags=10)
lag = 1
for lb, p in zip(lbvalues, pvalues):
    print(lag, lb, p)
    lag += 1

# 出力
1 4.095089120802025 0.04300796255246615
2 4.223295881654548 0.1210383379718042
3 6.047027807671336 0.10934450247698609
4 6.312955422660912 0.1769638685520115
5 6.457291061922424 0.26422887039323306
6 9.36446186985462 0.15409458108816818
7 40.47102807634717 1.0226763047711032e-06
8 45.2406378234939 3.312995830103468e-07
9 45.24892593869829 8.298135787344997e-07
10 46.35530787049628 1.2365386593885822e-06

lag が7未満の時は、p値が0.05を超えているので、帰無仮説$H_0$を棄却できず、
7以上の時は棄却できていることがわかります。

念の為ですが、lag が 8,9,10の時に棄却できているのは、
どれもデータが7点周期を持っていることが理由であり、
8,9,10点の周期性を持っていること意味するものではありません。

The Zen of Python

Pythonプログラマーならどこかでみたことがあると思われる「The Zen of Python」が、
pep 20 で定義されたものであることを知りました。

PEP 20 — The Zen of Python

ちなみに、読みたくなったら import thisでいつでも読むことができます。
thisって何だろうとずっと思っていたのですが、Easter Eggとして作られたものだったようです。


>>> import this
The Zen of Python, by Tim Peters

Beautiful is better than ugly.
Explicit is better than implicit.
Simple is better than complex.
Complex is better than complicated.
Flat is better than nested.
Sparse is better than dense.
Readability counts.
Special cases aren't special enough to break the rules.
Although practicality beats purity.
Errors should never pass silently.
Unless explicitly silenced.
In the face of ambiguity, refuse the temptation to guess.
There should be one-- and preferably only one --obvious way to do it.
Although that way may not be obvious at first unless you're Dutch.
Now is better than never.
Although never is often better than *right* now.
If the implementation is hard to explain, it's a bad idea.
If the implementation is easy to explain, it may be a good idea.
Namespaces are one honking great idea -- let's do more of those!

とても良いことが書かれていると思いますので、
定期的に読み返したいですね。

hyperoptで整数値の探索

今回もhyperoptの話題です。
この前の記事では、2変数関数の最小値を探すとき、
それぞれの変数を、一定区間の一様分布の中から探索しました。
そのときは、hp.uniformを使いました。

ただ、場合によっては、変数は実数ではなく整数を取ることもあります。
多くの例があると思いますが、機械学習で言えばDNNの中間層のユニット数などがそうですね。

この場合、 hyperoptで探索する空間を定義するときは、hp.randint(label, upper)というのを使うことで実現できます。
ドキュメントはこちら。
2.1 Parameter Expressions
ただ、これよく見ると最小値の指定ができません。
(ドキュメントに書かれてないだけで実は動くのではと思い試しましたが、本当にできません)

そのような時にどうしようかというのが今回の記事です。
0からではなく、1以上の数で探したいケースや、負の数も扱いたい時などありますからね。

今回使う関数はこちら。
$$f(x, y) = (x-5)^2 + (y+7)^2$$
この関数の$-10 \leq x\leq 10, -10 \leq y\leq 10$ の範囲で最小値を取る点を探します。
正解は明らかに$(x, y) = (5, -7)$ です。

それを実行するコードがこちら。


from hyperopt import hp, fmin, tpe, Trials


# 最小値を探したい関数
def f(x, y):
    return (x-5)**2 + (y+7)**2


# 探索する関数を定義する
def objective(args):
    x = args["x"]
    y = args["y"]

    return f(x, y)


# 探索する空間を定義する
space = {
        'x': -10 + hp.randint('x_', 21),
        'y': -10 + hp.randint('y_', 21),
    }

trials = Trials()

best = fmin(
            fn=objective,
            space=space,
            algo=tpe.suggest,
            max_evals=100,
            trials=trials,
        )

print(best)

# 以下出力
100%|██████████| 100/100 [00:00<00:00, 427.25it/s, best loss: 2.0]
{'x_': 16, 'y_': 2}

ちょっと惜しいですが、 f(x, y)が 2になる点を見つけたようです。

ポイントはここ


space = {
        'x': -10 + hp.randint('x_', 21),
        'y': -10 + hp.randint('y_', 21),
    }

このように定義することで、x_y_ はそれぞれ 0〜20の値を取り、
-10することで、xyは-10から10の値を取ります。
21にはならないので注意。(trials.vals["x_"]の値を確認するとx_=21は探していないことがわかります。)

出力の、{'x_': 16, 'y_': 2}から、
$x=16-10=6, y=2-10=8$ を見つけていることがわかります。

ただ、本当はx,yの値が欲しいのに、x_,y_を戻されると少しだけ不便です。
このようなときは、space_evalという関数が使えるようです。
ドキュメントが見つからないのですが、このページのサンプルコードで使われているのを参考にやってみます。


from hyperopt import space_eval
print(space_eval(space, best))

# 出力
{'x': 6, 'y': -8}

これで必要だった(x, y)の値が取れました。
正解から1ずつずれているのがちょっと惜しいですね。

hyperoptで多変数関数の最適化

前回前々回とhyperoptを扱ってきましたが、例があまりにもシンプルだったのでもう少しまともなサンプルを出します。
今回はあらかじめ定義された多変数関数が最大値をとる点を探します。
たった2変数の関数ですが、これ以上変数を増やしても本質的にはやることは変わりません。

きちんと最適解の近くを探していることがわかるように最後に可視化も行います。

例にする関数はこちらです。
$$2*e^{-(x-1)^2-(y-2)^2}+e^{-(x+3)^2-(y+4)^2}$$
点$(1, 2)$からちょっとずれたところで最大値を取りそうな関数ですね。
そして、 点$(-3, -4)$付近にもう一つ山があります。

早速やってみましょう。


import numpy as np
from hyperopt import hp,  fmin, tpe, Trials


# 最大値を探したい関数
def f(x, y):
    return 2 * np.exp(-(x-1)**2 - (y-2)**2) + np.exp(-(x+3)**2 - (y+4)**2)


# 探索する関数を定義する
def objective(args):
    x = args["x"]
    y = args["y"]

    # 最大値を検索するため、マイナスをつけて符号を反転
    return - f(x, y)


# 探索する空間を定義する
space = {
        'x': hp.uniform('x', -5, 5),
        'y': hp.uniform('y', -5, 5)
    }

trials = Trials()

best = fmin(
            fn=objective,
            space=space,
            algo=tpe.suggest,
            max_evals=100,
            trials=trials,
        )

print(best)

# 以下出力
100%|██████████| 100/100 [00:00<00:00, 293.11it/s, best loss: -1.9782848207515258]
{'x': 1.0647780122167085, 'y': 1.9180196819617374}

きちんとそれっぽい点を見つけてきましたね。
これまでの例と違うのは、最小値を探索する関数と、探索する空間が、fminの外側で事前に定義されている点です。

あとは trials オブジェクトの中身を確認し、正解付近を重点的に探索していることを確認しましょう。


import matplotlib.pyplot as plt
X = np.linspace(-5, 5, 100)
Y = np.linspace(-5, 5, 100)
xx, yy = np.meshgrid(X, Y)

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1)
ax.contour(xx, yy, f(xx, yy))
ax.scatter(trials.vals["x"], trials.vals["y"])
plt.show()

できた図がこちらです。

$(-3,-4)$付近はそこそこに、$(1, 2)$付近を重点的に探していることがわかります。

hyperoptで探索したパラメーターの履歴を確認する

前回の記事に続いて、hyperoptの記事です。
前回のサンプルコードでは、100回の探索を経て見つけた最良の値を取得しましたが、
そのほかの結果もみたいということはあると思います。
機械学習のモデルであれば、どのパラメーターが精度に寄与していたとか、
どれはあまり重要でないとか確認したいですからね。

そのような時、hyperopt では Trials というオブジェクトを使います。
ドキュメントのこちらを参考にやってみましょう。
1.3 The Trials Object


from hyperopt import fmin, tpe, hp, Trials
trials = Trials()
best = fmin(
            fn=lambda x: x ** 2,
            space=hp.uniform('x', -10, 10),
            algo=tpe.suggest,
            max_evals=100,
            trials=trials,
        )
print(best)

# 出力
100%|██████████| 100/100 [00:00<00:00, 554.94it/s, best loss: 8.318492974446729e-09]
{'x': 9.120577270352315e-05}

Trials() ってのが途中で出てきた以外は前回の記事と同じコードですね。
この時、trials変数に100回分の履歴が残っています。
試しに一つ見てみましょう。


# trials.trialsに配列型で格納されている結果の一つ目
trials.trials[0]

# 内容
{'state': 2,
 'tid': 0,
 'spec': None,
 'result': {'loss': 96.69903496469523, 'status': 'ok'},
 'misc': {'tid': 0,
  'cmd': ('domain_attachment', 'FMinIter_Domain'),
  'workdir': None,
  'idxs': {'x': [0]},
  'vals': {'x': [-9.833566746847007]}},
 'exp_key': None,
 'owner': None,
 'version': 0,
 'book_time': datetime.datetime(2019, 3, 10, 15, 4, 45, 779000),
 'refresh_time': datetime.datetime(2019, 3, 10, 15, 4, 45, 779000)}

$x=-9.833566746847007$ を試したら 損失関数$x^2$の値は$96.69903496469523$ だったということがわかります。

このほか、次のような変数にそれぞれの値の辞書も格納されています。
trials.results
trials.vals

せっかくなので、本当に正解の$x=0$付近を探索していたのか可視化してみましょう。


import matplotlib.pyplot as plt
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1)
ax.scatter(trials.vals["x"], trials.losses())
plt.show()

結果がこちら。

きちんと正解付近に点が密集していますね。

hyperoptのインストールと最も簡単な例

昨年参加したpycon2018
実践・競馬データサイエンス というセッションで知った、hyperopt というライブラリのメモです。
とりあえずこの記事では最もシンプルな例で動かすところまで行きます。
ちなみに、pyconの時の資料はこちら。
実践・競馬データサイエンス

なんでも、ハイパーパラメーターをグリッドサーチなどより効率的に探索してくれるとのこと。

公式サイトの通り、 pipでインストールできます。


pip install hyperopt

早速動かしてみましょう。
Basic tutorial

$x^2$ が最小になる$x$(正解は$x=0$です)を探索するコードです。


from hyperopt import fmin, tpe, hp
best = fmin(
            fn=lambda x: x ** 2,
            space=hp.uniform('x', -10, 10),
            algo=tpe.suggest,
            max_evals=100
        )
print(best)

# 以下出力
100%|██████████| 100/100 [00:00<00:00, 570.89it/s, best loss: 0.0009715421527694116]
{'x': 0.031169570942979175}

これで、 $-10<=x<=10$ の範囲を100回探索し、$x=0.031$あたりで最小と見つけてきたようです。 単純すぎていまいちありがたみが感じられませんね。 ただグリッドサーチと比較すると、仮に-10から10を単純に100分割したら区間の幅は0.2になるわけで、 それよりはずっと小さい誤差で最小値付近の値を見つけています。 (とはいえ、グリッドサーチなら $x=0$をピタリと見つけるとは思いますが。) 機械学習のハイパーパラメーターの探索等で役に立つと聞いて調べたものであり、 僕もそのように使うことが多いので、今後の記事でもう少し細かく紹介しようと思います。 一点気になるのは公式ドキュメントのこの記述。

Hyperopt has been designed to accommodate Bayesian optimization algorithms based on Gaussian processes and regression trees, but these are not currently implemented.

[翻訳]
Hyperoptは、ガウス過程と回帰木に基づくベイズ最適化アルゴリズムに対応するように設計されていますが、これらは現在実装されていません。

Qiitaなどで、hpyeroptをベイズ最適化のライブラリだって紹介している人を見かけたこともあるのですが、そうではなさそうです。

DataFrameを特定の列の値によって分割する

python(pandas)で巨大なデータフレームを扱っている時、
ある列の値によって、分割したいことが良くあります。

pandasのgroupbyのユーザーガイドを見ていて、良い方法を見つけたのでそれを紹介します。
Iterating through groups
のところに書かれていますが、groupbyした後そのままsumなどの関数を適用するだけでなく、イテレーターとして使うことで、
各グループに対して順番に処理を行えます。

実際にやってみましょう。
とりあえずサンプルのデータフレームを作ります。


import pandas as pd
import numpy as np
df = pd.DataFrame(
            np.random.randn(10000, 10),
            columns=['x'+str(i) for i in range(10)],
        )
df['y'] = np.random.randint(20, size=10000)

これで、長さが1万行、列が11個あるデータフレーム(df)ができました。
これをy列の値(0〜19)によって、20個のデータフレームに分割します。


df_dict = {}
for name, group in df.groupby('y'):
    df_dict[name] = group

これで、20個のデータフレームに分割できました。

ちなみに、この方法を知る前は次のように書いていました。


df_dict = {}
for name in set(df.y):
    df_dict[name] = df[df.y == name]

コードの量はほとんど同じなのですが、この方法は無駄な比較を何度も行うことになります。

実際、実行時刻を測ってみるとgroupbyを使う方法は、
CPU times: user 7.47 ms, sys: 3.28 ms, total: 10.8 ms
Wall time: 7.71 ms

なのに対して、使わない方法は
CPU times: user 16.4 ms, sys: 1.81 ms, total: 18.2 ms
Wall time: 16.5 ms
となり、倍以上の時間を要しています。

matplotlibで折れ線グラフ(plot)を書く時に線の書式を指定する

matplotlibで折れ線グラフを書く時に線の書式を指定する方法のメモです。
大抵のケースでは色で見分ければ十分なのですが、たまに点線や破線を使いたいことがあります。
(色が使えない場合もマーカーで見分けてもいいのですけどね。)

ドキュメントはこちらです。
set_linestyle(ls)
matplotlib.pyplot.plot の linestyle のリンクをクリックすると上のページに飛びます。

方法は簡単で、 plot する時に linestyle変数に ‘dashed’ などの文字列や ‘:’ などの記号を指定するだけです。
Noneで、線を消すこともできます。(その場合はマーカーを使いましょう)

試しに4種類の線を引いてみました。
(線の種類の区別を目立たさせるため、色は揃えました。)


import matplotlib.pyplot as plt
import numpy as np

# データ生成
x = np.arange(10)
y0 = np.random.normal(0, 0.5, 10)
y1 = np.random.normal(1, 0.5, 10)
y2 = np.random.normal(2, 0.5, 10)
y3 = np.random.normal(3, 0.5, 10)

plt.rcParams["font.size"] = 14
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(1, 1, 1, title="linestyleの例")
ax.plot(x, y0, linestyle="-", c="b", label="linestyle : solid, '-'")
ax.plot(x, y1, linestyle="--", c="b", label="linestyle : dashed, '--'")
ax.plot(x, y2, linestyle="-.", c="b", label="linestyle : dashdot, '-.''")
ax.plot(x, y3, linestyle=":", c="b", label="linestyle : dotted, ':''")
ax.legend(bbox_to_anchor=(1.05, 1))
plt.show()

出力はこちら。

pythonのdictから値を取得する時にデフォルト値を設定する

今は当たり前のように使っているテクニックですが、初めて知った時は感動したので紹介。

ドキュメントはここ

pythonの辞書オブジェクトから値を取得する時、[ ]を使う方法と、
get()を使う方法があります。


sample_dict = {
    'apple': 'リンゴ',
    'orange': 'オレンジ'
}

print(sample_dict['apple'])  # 'リンゴ'
print(sample_dict.get('orange'))  # 'オレンジ'

タイプ数が少ないので[ ]の方を好むかたが多いと思うのですが、
getの方には、引数がkeyの中に存在しなかった時の初期値を指定できるメリットがあります。

例えば、
sample_dict["grape"] は keyErrorが発生してしまいます。
しかし、
sample_dict.get("grape")はNoneを返しますし、
sample_dict.get("grape", "未定義")とすると、
2個目の引数である、”未定義”を返してくれます。

エラー処理もいらないし、事前に辞書に取得しようとしているkeyが含まれているか確認する必要もなくなり、
非常に便利です。

requestsにタイムアウトを設定する

requestsを使ってurlのリストをチェックしていた時にredirectに加えて困ったのが、応答に時間がかかるサーバーへアクセスした時です。
数秒待って結果が戻って来ればまだいいのですが、そのままスクリプトが進まなくなってしまうことがありました。
これを防ぐには、timeoutを設定すると良いようです。
(デフォルトでは設定されていない)

ドキュメントはこちら。
クイックスタート – timeouts
日本語が少しおかしい気がします。
timeout パラメーターに秒数を指定すると、指定した秒数の間、Requestsのレスポンスの待機を止めることができます。
おそらく意図は「timeout パラメーターに秒数を指定すると、指定した秒数でRequestsのレスポンスの待機を止めることができます。」ではないかと。

早速ですがやってみましょう。
timeout になると例外が発生するので、キャッチできるようにします。
(本来はrequestsライブラリで定義されている専用の例外を使うべきなのですが、とりあえずException使います)


import requests

url = "https://httpstat.us/200?sleep=10000"  # 時間のかかるURL

try:
    response = requests.get(url, timeout=3)
    print(response.status_code)  # 実行されない
    print(response.url)  # 実行されない

except Exception as e:
    print(e.args)

# 出力
(ReadTimeoutError("HTTPSConnectionPool(host='httpstat.us', port=443): Read timed out. (read timeout=3)",),)

想定通りに動きました。