定常過程の例としてのiid過程とホワイトノイズ

前回の記事で定義した定常過程の例を紹介します。
(今回も沖本本からの引用です。)

iid系列

iid系列は、強定常過程の例です。

定義:
各時点のデータが互いに独立でかつ同一の分布に従う系列はiid系列と呼ばれる。

iidとは、 independent and identically distributedの略です。
Wikipedia : 同時独立分布

$y_t$が期待値$\mu$,分散$\sigma^2$のiid系列であることを $y_t\sim iid(\mu,\sigma^2)$と書きます。

ホワイトノイズ

独立性や同一分布性は非常に強い仮定なのですが、
実用上、モデルの錯乱項などとして使うのであればもう少し弱い仮定で十分です。
そこで次のホワイトノイズというものが使われます。
ホワイトノイズは弱定常過程です。

定義:
すべての時点$t$において
$$
\begin{align}
E(\varepsilon _t) &=0\\
\gamma _k &= E(\varepsilon _t \varepsilon _{t-k}) =\left\{
\begin{array}{ll}
\sigma^2, & (k=0)\\
0, & (k\neq 0)
\end{array}\right.
\end{align}
$$
が成立するとき, $\varepsilon _t$はホワイトノイズ(white noise)と呼ばれる.

$\varepsilon _t$が分散$\sigma^2$のホワイトノイズであることを $\varepsilon _t\sim W.N.(\sigma^2)$と書きます。
次のようにホワイトノイズを使って、基礎的な弱定常過程の例を作れます。
$$
y_t = \mu + \varepsilon _t,\ \ \ \varepsilon _t\sim W.N.(\sigma^2)
$$

自分だけかもしれないのですが、ホワイトノイズと聞くと無意識に$(\varepsilon _t,\cdots,\varepsilon _{t+k})$の同時分布が期待値0の多変量正規分布であるかのようにイメージしてしまいます。
しかし定義の通り、ホワイトノイズにはそこまでの強い仮定はないので注意です。

ちなみに $(y_t,\cdots,y_{t+k})$の同時分布が多変量正規分布の場合、それは正規過程と呼ばれ、
正規過程に従うホワイトノイズは正規ホワイトノイズと呼ばれます。(沖本本のp9,p12)
正規ホワイトノイズはiid過程なので、強定常過程になります。

確率過程の弱定常性と強定常性

時系列データの分析において重要な定常性の定義のメモです。
用語は沖本先生の経済・ファイナンスデータの計量時系列分析の定義に従っています。

確率過程の定常性とは、その過程の同時分布や基本統計量の時間不変性に関するものです。
何を不変とするかによって、弱定常性 (weak stationarity)と強定常性(strict stationarity)に分類されます。

以下、$y_t$を確率過程とします。

弱定常性

弱定常性は、過程の期待値と自己共分散が一定であることです。

定義:
任意の$t$と$k$に対して、
$$
\begin{align}
E(y_t) & = \mu\\
Cov(y_t, y_{t-k}) & = E[(y_t-\mu)(y_{t-k}-\mu)] = \gamma_k
\end{align}
$$
が成立する場合、過程は弱定常(weak stationary)と言われます。

$y_t$の期待値は任意の$t$に対して一定で、$y_t$と$y_{t-k}$の共分散は$k$の値によってのみ決まります。
また、この二つが時点に依存しないので、自己相関も時点に依存しません。

経済・ファイナンスの分野では単に定常性というと、弱定常性をさすことが多いそうです。(沖本本 P10)

強定常性

強定常性は、同時分布が不変であることを要求します。

定義:
任意の$t$と$k$に対して、$(y_t,y_{t+1},\cdots,y_{t+k})$の同時分布が同一となる場合、
過程は強定常(strict stationary)と言われます。

弱定常は期待値と共分散だけに対する仮定だったので、強定常のほうが名前の通り強い仮定になります。

ただし、強定常ならば絶対に弱定常というわけではなく、
“過程の分散が有限であるならば、”強定常過程は弱定常過程である、
という風に一つ条件が必要です。

pandasで時系列データの自己相関係数算出

時系列データとそれをn時間分シフトさせたデータ間の相関係数を自己相関係数と呼びます。
データに周期性があるかどうかを確認するときなどに調査します。

単に配列をスライドさせてnumpyか何かで相関係数を出してもいいのですが、
pandasに専用の関数が定義されているのでそれを使って求めることもできます。

pandas.Series.autocorr

試しにやってみましょう。


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

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

# lag=0から29までの自己相関係数
auto_series = [series.autocorr(lag=i) for i in range(30)]

# 可視化
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(2, 1, 1, title="周期性のあるデータ")
ax.plot(series.index, series)
ax = fig.add_subplot(2, 1, 2, title="自己相関係数")
ax.bar(range(30), auto_series)
plt.show()

出力されるグラフがこちら。

7点間隔の周期が綺麗に現れているのがわかります。

LaTeXで行列の成分を省略したときのドットを書く

大学院生のときは頻繁に書いていたのに、すっかり忘れてしまっていたので$\LaTeX$の復習です。

まず、要素を省略せずに3*3行列を書く例。

$$
\textbf{A} =
\left(
\begin{array}{ccc}
a_{11} & a_{12} & a_{13} \\
a_{21} & a_{22} & a_{23} \\
a_{31} & a_{32} & a_{33}
\end{array}\right)
$$

結果。
$$
\textbf{A} =
\left(
\begin{array}{ccc}
a_{11} & a_{12} & a_{13} \\
a_{21} & a_{22} & a_{23} \\
a_{31} & a_{32} & a_{33}
\end{array}\right)
$$

次に、ドットを使って書く例。
横並びは\cdots, 縦並びは\vdots, 斜めは\ddots を使います。


$$
\textbf{W} =
\left(
\begin{array}{ccc}
w_{11} & \cdots & w_{1n} \\
\vdots & \ddots & \vdots \\
w_{n1} & \cdots & w_{nn}
\end{array}\right)
$$

結果。
$$
\textbf{W} =
\left(
\begin{array}{ccc}
w_{11} & \cdots & w_{1n} \\
\vdots & \ddots & \vdots \\
w_{n1} & \cdots & w_{nn}
\end{array}\right)
$$

scipyで手軽にカイ二乗検定

業務でカイ二乗検定を行う場面は、割と多くあります。
自分の場合は、ABテストの評価で使うことが多いです。

Excelでも自分でコードを書いてでもさっとできるのですが、
scipyに専用の関数が用意されているのでその紹介です。

scipy.stats.chi2_contingency

これを使うと、引数にデータのクロス表を放り込むだけで、
カイ二乗値、p値、自由度、期待度数表を一発で算出してくれます。
3つの値と、期待度数がarray形式ではいったタプルが戻るので、引数を4つ用意して受け取ると便利です。


import pandas as pd
from scipy.stats import chi2_contingency

# 検定対象のデータ(サンプルなので値はダミー)
df = pd.DataFrame([[300, 100], [800, 200]])

chi2, p, dof, expected = chi2_contingency(df, correction=False)

print("カイ二乗値:", chi2)
print("p値:", p)
print("自由度:", dof)
print("期待度数:", expected)

# 以下出力
カイ二乗値: 4.242424242424244
p値: 0.03942583262944296
自由度: 1
期待度数: [[314.28571429  85.71428571]
 [785.71428571 214.28571429]]

このようにめんどな手間がなく一瞬で検定ができました。
ちなみに、
correction=False
はイエーツの修正を”行わない”設定です。
デフォルトがTrueなので注意してください。

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ずつずれているのがちょっと惜しいですね。

Data Driven Developer Meetup #5 に参加しました

先日、Data Driven Developer Meetup #5 に参加してきました。
#4にも参加しているので2回目です。

今回の会場はログリーさんのオフィス。職場から近かったので徒歩で行けて便利でした。

メインセッション

データドリブンな小売店舗経営を支えるプロダクト開発の裏側』 株式会社ABEJA 大田黒 (@xecus) さん

まずこの発表資料は ABEJA SIX 2019で使われたものだそうです。
ABEJA SIXは参加したかったのですが業務的に行けなかったので、ラッキーでした。

動画の顔認識や、店内の導線の可視化など技術力が凄まじく高いのはもちろんですが、
「リアル版のグーグルアナリティクス」といった非常にイメージをつかみやすい説明が印象的でした。
これはビジネス的に強そうです。
上のスライド資料だと動画が静止画になってしまいっていましたが、
店舗内入り口のカメラの動画で、性別と年齢が推定されているデモなどは驚きでした。

プロジェクト進行における問題のところは自分も似たような問題に直面したのでよくわかります。
ゴール設定の難しさも、不確実性の多さも、本当にウォーターフォールとは相性が悪いです。

まだ自分の職場では、「研究開発」と言えるほどモデルの開発に専念できる体制も作れておらず、
機械学習のプロジェクトはもっとごちゃごちゃしているのですが、
チームで進める土壌づくりは力を入れないといけないと思いました。

機械学習で動くシステムを実運用に組み込むということ』 株式会社ブレインパッド アナリティクスサービス部 吉田勇太 (@yutatatatata ) さん

これもあるあるのオンパレードで苦笑いする場面の多い発表でした。
画像とテキストという違いはありますが、目視でやっていたことを機械学習で行いたいというニーズは自分も経験しています。
その対象が不均衡データで、
PrecisionとRecallがトレードオフになる厄介さや、Recallが重要だからと言って、Precisionはどこまでも犠牲にしていいわけじゃない苦労など
経験してきたので非常によくわかりました。

また、良いモデルを作るだけではなく業務フローをしっかり考えないといけないこと、
それにあたって、他のメンバーの協力を得ることの重要性など自分の失敗経験もあるので本当に納得です。

このあと、LTはお酒も入って気楽な雰囲気で聞いていました。

LT

学習行動データを蓄積するLearnig Record Storeの開発事例』by @yukinagae さん

まだまだ取り組みが始まったばかりという印象が大きかったのですが、
「Learnig Record Store」というコンセプトをしっかり立てて進められているのが素晴らしいと思いました。
GCPUG覚えておきます。

データドリブンを提供するサービスBrand Officialのアーキテクチャについて』 by @YasutakaYamamoto

久々に Airflowではなくdigdagの方が登場しました。
まだdigdagは使ったことがないのですが一度くらいは調べておきたいですね。
embulkは使っているのですが、よく理解していないので調べなければ。

今回も非常に参考になること、励みになることが多く参加してよかったです。
運営の皆様、登壇者の皆様、ありがとうございました。

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()

結果がこちら。

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