ベクトル自己回帰モデル

久しぶりに時系列分析の話題です。
今回はベクトル自己回帰(VAR)モデルを紹介します。
参考書はいつもの沖本先生の経済・ファイナンスデータの計量時系列分析です。
(VARモデルは第4章)

ベクトル自己回帰(VAR)モデルは、自己回帰(AR)モデルをベクトルに一般化したものです。
自然数$p$に対して、ベクトル $\mathbf{y}_t$ を定数(のベクトル)と、自身の$p$期以内の過去の値に回帰したモデルになります。
つまり、次のようなモデルになります。
$$
\mathbf{y}_t = \mathbf{c}+\mathbf{\Phi}_1\mathbf{y}_{t-1}+\cdots+\mathbf{\Phi}_p\mathbf{y}_{t-p}+\mathbf{\varepsilon}_t,
\ \ \ \mathbf{\varepsilon}_t\sim W.N.(\mathbf{\Sigma})
$$
($\mathbf{\varepsilon}$がなぜか太字にならない。自分の環境だけかな。)

ここで、$\mathbf{c}$ は $n\times 1$定数ベクトルで、$\mathbf{\Phi}_i$は$n\times n$行列です。
$\mathbf{\varepsilon}_t$はベクトルホワイトノイズと呼ばれるものです。
(単純に各要素がホワイトノイズのベクトルというわけではありません。このブログでは新出語なので後日別の記事で取り上げます。)

2変量のVAR(1)モデル、要するに実質的なVARモデルとしては一番小さいモデルを具体的に表記すると次のようになります。

$$\left\{\begin{matrix}\
y_{1t} &=&c_1+\phi_{11}y_{1,t-1}+\phi_{12}y_{2,t-1}+\varepsilon_{1t},\\\
y_{2t} &=&c_2+\phi_{21}y_{1,t-1}+\phi_{22}y_{2,t-1}+\varepsilon_{2t}\
\end{matrix}\
\right.\
\left(\begin{matrix}\varepsilon_{1t}\\\varepsilon_{2t}\end{matrix}\right)\sim W.N.(\mathbf{\Sigma})$$

$$
\Sigma = \left(\begin{matrix}\
\sigma_1^2 & \rho\sigma_1\sigma_2\\\
\rho\sigma_1\sigma_2 & \sigma_2^2\
\end{matrix}\right)\
$$

さて、具体的にモデルの形を書き出したのでパラメーターの数を数えてみたいと思います。
回帰式の係数(4個)と定数(2個)で計6個、さらに分散共分散行列部分が3個で、合計9個のパラメーターを持っていることがわかります。

一般的に$n$変数のVAR(p)モデルの場合、
回帰式の定数が$n$個、係数が$n^2p$個で合計$n(np+1)$個、分散共分散行列部分が$n(n+1)/2$個のパラメーターを持ちます。
そのため、比較的大きなモデルになってしまいます。

自分の経験でも、モデルに組み込みたい変数の値の数($n$)は10個も20個もあって、ラグ($p$)は3期くらい前まで見たい、
ということがありましたが、集められたデーターが全然足りず、推定が全然できないということがありました。
(そのときは結局変数をかなり絞り込んだモデルをいくつも作って個別に検証しました。)

さて、ベクトル自己回帰モデルを使う目的ですが、主に二つあります。
一つは、予測の精度向上です。
自己回帰モデルでは、ある値が自身の過去の値で回帰できることこを仮定して予測しますが、
ある値が自分の過去の値以外の外部要因に一切影響を受けない、ということは実用上あまりなく、
大抵は複数の値が相互に関係しているので、複数の種類の値を用いて予測することにより精度の向上が期待できます。

もう一つは、複数の変数間の関係を分析することです。
モデルの形からも明らかなように、それぞれの変数が他の変数にどれだけ寄与しているかをみることができます。

pandasで縦横変換(pivot_table)

前回の更新でPrestoでデータの縦横変換をする方法を紹介しましたが、
クエリで処理を完結させる必要がないときは、一旦pandasのデータフレームに格納してから処理をするのも便利です。

その場合、 pandas.pivot_table を使います。

使い方は簡単で、pd.pivot_tableに、変換したいデータフレーム、
列にするカラム、行にするカラム、集計する値、集計に使う関数を指定するだけです。
fill_value引数で欠損値を埋めるなどの細かい設定もできます。
ドキュメントの例を使ってやってみます。


import pandas as pd
# データ作成
df = pd.DataFrame(
    {
        "A": ["foo", "foo", "foo", "foo", "foo", "bar", "bar", "bar", "bar"],
        "B": ["one", "one", "one", "two", "two", "one", "one", "two", "two"],
        "C": ["small", "large", "large", "small",
              "small", "large", "small", "small", "large"],
        "D": [1, 2, 2, 3, 3, 4, 5, 6, 7],
        "E": [2, 4, 5, 5, 6, 6, 8, 9, 9]
    }
)
print(df)
"""
     A    B      C  D  E
0  foo  one  small  1  2
1  foo  one  large  2  4
2  foo  one  large  2  5
3  foo  two  small  3  5
4  foo  two  small  3  6
5  bar  one  large  4  6
6  bar  one  small  5  8
7  bar  two  small  6  9
8  bar  two  large  7  9
"""

table_0 = pd.pivot_table(
                df,
                values="D",
                index="A",
                columns="C",
                aggfunc="sum",
                fill_value=0,
        )
print(table_0)
"""
C    large  small
A
bar     11     11
foo      4      7
"""

# 行や列、集計関数は配列で複数指定することもできる
table_1 = pd.pivot_table(
                df,
                values="D",
                index=["A", "B"],
                columns="C",
                aggfunc=["sum", "count"],
                fill_value=0,
        )
print(table_1)
"""
          sum       count
C       large small large small
A   B
bar one     4     5     1     1
    two     7     6     1     1
foo one     4     1     2     1
    two     0     6     0     2
"""

PrestoのMap型を使った縦横変換

前回の記事で紹介したPrestoのMap型ですが、これを使うとデータの縦横変換(ピボット)がスマートに行えます。

参考: ピボットテーブル&チャート

上記リンク先のトレジャーデータの記事中の画像のテーブルを例に説明します。

縦持ちのテーブル (vtable)
uid key value
101 c1 11
101 c2 12
101 c3 13
102 c1 21
102 c2 22
102 c3 23

このようなテーブルを次のように変換したいとします。

横持ちのテーブル
uid c1 c2 c3
101 11 12 13
102 21 22 23

この変換は、MAP_AGGを使って、key列とvalue列の対応のMapを一旦作成し、
それぞれのMapから各Keyの値を取り出すことで実現できます。
具体的にクエリにしたのが次です。


SELECT
    uid,
    kv['c1'] AS c1,
    kv['c2'] AS c2,
    kv['c3'] AS c3
FROM (
    SELECT
        uid,
        MAP_AGG(key, value) AS kv
    FROM
        vtable
    GROUP BY
        uid
) AS t

個人的には服問い合わせはあまり好きではなく、PrestoではWITHを使って書きたいので次のようにすることが多いです。


WITH
    t AS (
        SELECT
            uid,
            MAP_AGG(key, value) AS kv
        FROM
            vtable
        GROUP BY
            uid
    )
SELECT
    uid,
    kv['c1'] AS c1,
    kv['c2'] AS c2,
    kv['c3'] AS c3
FROM
    t

Prestoではない、(MAP型のない)通常のSQLでは次のように書かないといけないのですが、
上記のMap型を使ったものの方が随分すっきりかけているように見えます。
(MAP_AGGを知らない人には読めないのが難点ですが)


SELECT
    uid,
    MAX(
        CASE WHEN key = 'c1' THEN
            value
        ELSE
            NULL
        END
    ) AS c1,
    MAX(
        CASE WHEN key = 'c2' THEN
            value
        ELSE
            NULL
        END
    ) AS c2,
    MAX(
        CASE WHEN key = 'c3' THEN
            value
        ELSE
            NULL
        END
    ) AS c3
FROM
    vtable
GROUP BY
    uid

PrestoのMap型について

Prestoには Mapというデータ型があります。
これは Pythonのdictと似たようなデータ型で、 キーと値のペアで構成されるものです。

ドキュメントを見ると、
MAP(ARRAY['foo', 'bar'], ARRAY[1, 2]) というMAPを作る関数の例が紹介されていますが、実行すると次のような結果を得ることができます。


-- クエリ
SELECT
    MAP(ARRAY['foo', 'bar'], ARRAY[1, 2])

-- 以下結果
{"bar":2,"foo":1}

以前紹介した、 TD_PARSE_AGENT という関数がありますが、この結果もMapです。
参考: TreasureDataのTD_PARSE_AGENT関数が便利

Map から、キーをしていて値を取得するときは、上のTD_PARSE_AGENTの記事でやっているように、map に [‘key’]をつけて取得するか、
ELEMENT_AT(map(K, V), key) を使います。
個人的にはPythonなどと近い書き方の [ ] を使う方法の方が好きです。

そのほかの、Mapの使い方ですが、ドキュメントの
6.17. Map Functions and Operators
というページにまとまっているのでこちらがわかりやすいです。

さて、 Tableのある列の値を キーとし、別の列の値を値とするMapを作りたくなる場合があります。
そのようなときは、 MAP_AGG という関数が使えます。
ドキュメントでは別のページにあるので探しにくいのですが、こちらにあります。
6.14. Aggregate Functions
の Map Aggregate Functions。

キーに指定した列に重複があったらどうなるかが気になったのですが、
試してみたところ、バリューの列のどれかの値が何らかの規則で一つだけ選ばれて採用されるようです。
(どのような基準で選ばれているのかは結局わかりませんでした。使うときは注意が必要ですね。)

キーに重複があり、バリューを(配列の形で)全部残したいときは MULTIMAP_AGG が使えます。

MAP_AGG で MAPを作って、 map[key] で値にアクセスする、ということだけ覚えておけば、
特に問題なく使うことができると思います。

2020年のご挨拶と今年の目標

新年明けましておめでとうございます。本年もよろしくお願いします。

年明け最初の投稿なので、本年の目標をまとめておきたいと思います。
細かい目標は他にも多々ありますが、データサイエンティストとしてのスキルアップの観点では、
次の二つを重点目標として、取り組みます。

1. kaggleのコンペに挑戦する

一つ目の目標はこれです。今年はブログ以外の技術アウトプットとして、kaggleに挑戦したいと思っています。
データサイエンティストではありますが、業務では機械学習以外の仕事が多く、
ぼーっとしているとなかなか技術が伸びていきません。
そのままだと稀に機械学習のタスクが発生した時に十分な成果が出せず非常に苦労します。
今は東京大学の通信講座であるDL4USを受講していますが、もうまもなく最終課題も終わるので、
その次の挑戦としてコンペに参加していこうと思います。
もちろんメダルを目指したいですが、まずはコンペに取り組む習慣づくりからはじめて、
常に何かしらのコンペに取り組んでいる状態を維持します。

2. ブログ記事100件投稿

昨年は306記事投稿していますが、予告通りペースを落とします。
ただ、それでも投稿自体は継続し、1週間に2記事のペースで年間100記事を目指します。
内容も見直していく予定です。
いきなりキャリア系のポエム記事ばかりにするつもりはなく、今まで通りのノリの投稿も続けますが、
キャリア系の記事でも誰かの参考になることがあると思うので、徐々に発信の幅を広げていきたいです。
特に30歳を超えて未経験の状態からいきなりデータサイエンティストに転職してきた頃の話とか、
その後、分析チームのマネジメントや採用業務等も担当するようになっているので、
何かしら誰かの参考になる話もあるかと思います。

2019年のまとめ

2019年最後の更新です。
今年から始めたブログでしたが無事に1年間更新を続けることができました。
ということでGoogleアナリティクスのデータ等も見ながら、振り返って見ます。
(記事執筆時点のデータなので後日残り1日文のデータを入れて更新するかも)

まず基本的なデータ。

– 記事数 305 記事 (この記事含む)
– 訪問ユーザー数 23,803人
– ページビュー 36,569回

密かに目標にしていた300記事は無事に達成し、アクセスもそこそこ集まるようになりました。
どこまで役に立っているのかわかりませんが少なくとも更新する意味はあるブログになってきたのかなと思います。

実際、更新してきてどうだったか、という観点の話はつい先日300記事達成の投稿で書いたので省略します。

次は今年よく読まれた記事の紹介です。
トップテンは次のようになりました。

  1. macにgraphvizをインストールする
  2. pythonでARモデルの推定
  3. DataFrameを特定の列の値によって分割する
  4. graphvizで決定木を可視化
  5. pandasでgroupbyした時に複数の集計関数を同時に適用する
  6. pythonで編集距離(レーベンシュタイン距離)を求める
  7. scikit-learnでテキストをBoWやtfidfに変換する時に一文字の単語も学習対象に含める
  8. Prestoで1ヶ月後の時刻を求める時に気をつけること
  9. pythonで累積和
  10. pythonでARMAモデルの推定

ベスト5は3Qの時とほとんど変化していませんが、6位から10位は少し意外なのも入っていますね。累積和など。
データサイエンティスト色をもっと強めたいなとは思うのですが、とりあえず技術ブログっぽいものにはなったと思います。

元々、各所に散らばってしまっていた自分のメモや検証結果などを一箇所にまとめたいという思いで始め、
訪問者のことよりも自分にとっての使いやすさ重視の記事が多いブログですが、
多くのかたに訪問していただきとてもありがたく思います。

さて、来年の更新ですが、今のペースでやっていくのは少し難しいと思っているので、更新頻度は見直したいと思っています。
というのも、最近は自分の学習時間においてブログへのアウトプットの比重が高まりすぎ、
腰を据えたインプットがおろそかになっているという課題も感じているからです。

年末年始に来年の目標と計画を整理し、このブログの運用はちょうどいい塩梅を探しながら続けていきたいと思います。
とりあえずお正月期間は更新をお休みし、来年の更新は6日以降から再開予定です。

今年一年ありがとうございました。良いお年を。

Googleアナリティクスでアラート設定

このブログでもアクセスの分析のためにGoogleアナリティクスを使っています。
順調に訪問者数が増えているのですが、何百人/日といった節目の達成はやはり嬉しいので、
数ヶ月前から、達成したら通知が来るように設定しています。

そのために使っているのがカスタムアラート機能です。

ドキュメント: カスタム アラートの作成、管理

左ペインのカスタムの中にあるカスタムアラートから飛べます。(飛び先は管理の中なのでそちらからも遷移できます)

アラート名を設定し、
期間を 日/週/月から選択、
このアラートが発生したときにメールで通知する。 にチェックを入れてメアドを設定、
あとは適用するトラフィックの対象と、指標、閾値となる数値を入れて保存したら完成です。

このブログだと、全てのトラフィックで、ユーザー数などを見ているだけですが、
ページを絞った直帰率や新規セッション割合など、かなり多様な設定が可能です。

個人ブログであればアラートが上がるとモチベーションが上がるようなものをいくつか設定しておくと良いですし、
仕事で分析しているサイトであれば異常時の通知などのに使えると思います。

関数内で発生した例外を呼び出し元にも伝える

昨日に続いて例外処理の話です。
ある関数内に、エラーが発生しうる処理がある時、その関数内でtry:〜except:〜処理を書いて
綺麗に例外を処理するけど、その例外を関数の呼び出し元にも伝えて例外を伝播させたいことがあります。

このような時も、 raise 文を使うことができます。

自分は最近まで raise 文は新規に例外を発生させいさせる機能しかないと思ってました。
こういう風に。


raise ValueError

"""
ValueError                                Traceback (most recent call last)
 in ()
----> 1 raise ValueError

ValueError: 
"""

しかし、raiseを単体で使用すると、そのスコープで有効になっている例外を再送出できます。

参考: 7.8. raise 文

試してみる前に、非常に単純な例なのですが次のようなケースを考えてみます。
関数 inv は 引数の逆数を返す関数で、0が渡されたら例外になるはずのものです。
そして、0が渡されたら内部で例外処理をしています。
そして、print_inv は inv を使って、与えられた数の逆数を表示します。


def inv(x):
    try:
        return 1/x
    except Exception as e:
        print(e)


def print_inv(x):
    try:
        print(inv(x))
    except Exception as e:
        print(e)
    else:
        print("print_invで例外は発生しませんでした")


print_inv(0)
"""
division by zero
None
print_invで例外は発生しませんでした
"""

ご覧の通り、 inv 内で例外処理しているので、呼び出し元では例外の発生を検知できていません。

ここで、 raiseを使って例外の再送出を入れてみます。


def inv(x):
    try:
        return 1/x
    except Exception as e:
        print(e)
        # 例外の再送出
        raise


def print_inv(x):
    try:
        print(inv(x))
    except Exception as e:
        print(e)
    else:
        print("print_invで例外は発生しませんでした")


print_inv(0)
"""
division by zero
division by zero
"""

division by zero が 2回表示されました。 inv と print_inv でそれぞれキャッチされた例外です。

Pythonにおける例外処理

jupyterでインタラクティブにPythonを使っているとあまり必要ないのですが、
本番コードを書くときなどは流石に例外処理を真面目に実装する必要があることがあります。
そこまで高頻度にあることではなく、すぐ忘れてしまうので、書き方をまとめておこうと思います。

参考になるドキュメントは次の2箇所です。
8. エラーと例外
組み込み例外

基本的に次のような書き方になります。
必須なのは、 try と except で、 exceptは複数書くこともできます。
except する例外には as e のように別名をつけることができ、
別名をつけておけば処理中で利用できます。
else と finally はオプションなので不要ならば省略可能です。


try:
    # ここに例外が発生しうるコードを書く

except [キャッチしたい例外クラス]:
    # 例外が発生した時に実行するコード

else:
    # 例外が発生なかった時に実行するコード

finally:
    # 必ず実行するコード

とりあえず定番の 0で割る演算で試してみましょう。


import numpy as np


def inv(data):
    try:
        inverse_data = 1/data
    except ZeroDivisionError as e:
        print(e)

    except TypeError as e:
        print(e)

    else:
        print("正常終了")
        return inverse_data
        print("このメッセージは表示されない")

    finally:
        print("finallyに書いた文は必ず実行されます")


print(inv(5))
"""
正常終了
finallyに書いた文は必ず実行されます
0.2
"""

print(inv(0))
"""
division by zero
finallyに書いた文は必ず実行されます
None
"""

print(inv("a"))
"""
unsupported operand type(s) for /: 'int' and 'str'
finallyに書いた文は必ず実行されます
None
"""

例外が発生した、0と”a” については想定通りに動きました。

実は例外が発生しなかったinv(5)が僕にとっては少し驚きでした。
else: のブロックの中で、 return して関数を抜けているので、
それより後ろの finally: のブロックは流石に実行されないと思っていたのですが、
print関数がバッチリ実行されています。

改めてよく読んでみれば、ドキュメント中にもしっかりそう書いてありました。
この辺りはきちんと理解して使う必要がありそうです。

– もし try 文が break 文、 continue 文または return 文のいずれかに達すると、その:keyword:break 文、 continue 文または return 文の実行の直前に finally 節が実行されます。
– もし finally 節が return 文を含む場合、 try 節の return 文より先に、そしてその代わりに、 finally 節の return 文が実行されます。

今回はブログ記事用に書いたコードだったので、
ZeroDivisionError と TypeError を 分けて書きましたが、
Exception のようなキャッチできる範囲の広い例外を指定しておけばまとめて受け取ってくれます。
(本当はあまり良くないと思うのですが、便利なので大抵そうしています。)


def inv2(data):
    try:
        inverse_data = 1/data
        return inverse_data
    except Exception as e:
        print(e)
        return None


print(inv2(5))
"""
0.2
"""

print(inv2(0))
"""
division by zero
None
"""

print(inv2("a"))
"""
unsupported operand type(s) for /: 'int' and 'str'
None
"""

また、例外処理の中で例外の種類を区別する必要が全くない場合、
except Exception as e:
の代わりに、
except:
とだけ書いておけば、より簡単に全ての例外をキャッチしてくれます。

DataFrameをマージする時にkeyの一意性を確認する

昨日の indicator の記事を書くためにドキュメントを読んでいて見つけた、 validate という引数の紹介です。

データフレーム通しを結合する時に、結合に使うキーのユニーク性が重要になることがあります。
事前に確認するようにコードを書いておけば済む話ではあるのですが、
pandasのmerge関数では、キーが一意でなかった時にエラーを上げてくれる引数があるようです。

ドキュメント: pandas.merge

引数 validate には、 one_to_one / one_to_many / many_to_one / many_to_many の4種類の文字列か、None(デフォルト)を
渡すことができます。
そして、 one_to_one なら 1:1, one_to_many なら 1:m 、という風にkeyが対応してなければエラーを上げてくれます。
many_to_many と None はノーチェックです。
なので、きちんとエラーをキャッチするように次のようにコードを書けます。
(データの準備等は省略)


try:
    df_merge = pd.merge(
            df_0,
            df_1,
            on="key",
            how="outer",
            validate="one_to_many",
        )
except pd.errors.MergeError as e:
    print(e)

# Merge keys are not unique in left dataset; not a one-to-many merge

発生するエラーは pd.errors.MergeError です。
データを何種類か用意して、validateの引数を変えながら動かすと色々動きがわかると思います。