Pythonで多変量正規分布に従う乱数を生成する

ベクトル自己回帰のダミーデータを生成するために、多変量正規分布に従う乱数が必要なので、
Pythonで生成する方法を紹介します。

numpyとscipyにそれぞれ用意されています。
同じ名前の関数だったので、どちらかの実装をもう一方がラップしているのかと思っていたのですが、
引数の微妙な違いなどあり、どうやら個別に実装されているようです。

ドキュメントはそれぞれ次のページにあります。

numpy
numpy.random.multivariate_normal
(この記事の numpy は version 1.16を使っています。 numpy 1.17.0 のリリースノートを見ると、random moduleに変更が加えられており、どうやらこの関数にも影響が出てるようなのでご注意ください。)

scipy
scipy.stats.multivariate_normal

さて、実際に期待値と分散共分散行列を指定してそれぞれ乱数を生成してみましょう。


import numpy as np
from scipy.stats import multivariate_normal

# 期待値と分散共分散行列の準備
mean = np.array([3, 5])
cov = np.array([[4, -1.2], [-1.2, 1]])

# numpy を用いた生成
data_1 = np.random.multivariate_normal(mean, cov, size=200)

# データ型の確認
print(data_1.shape)
# (200, 2)

# 期待値の確認
print(np.mean(data_1, axis=0))
# [3.00496708 4.94669956]

# 分散共分散の確認
print(np.cov(data_1, rowvar=False))
"""
[[ 3.86542859 -1.31389501]
 [-1.31389501  0.93002097]]
"""

# scipyで生成する方法
data_2 = multivariate_normal(mean, cov).rvs(size=200)

# データ型の確認
print(data_2.shape)
# (200, 2)

# 期待値の確認
print(np.mean(data_2, axis=0))
# [2.81459692 5.10444347]

# 分散共分散の確認
print(np.cov(data_2, rowvar=False))
"""
[[ 4.46151626 -1.28084696]
 [-1.28084696  1.06831954]]
"""

それぞれきちんと生成できたようです。

分散共分散行列の正定値性のバリデーションなど細かなオプションを持っていますが、
あまり使う機会はなさそうです。
(きちんと行う場合も、事前に固有値を求めて確認しておけば大丈夫だと思います。)

ベクトルホワイトノイズ

前回の記事で弱定常過程の定義をベクトルに拡張しましたが、今回はホワイトノイズをベクトルに拡張した
ベクトルホワイトノイズを定義します。
参考: 定常過程の例としてのiid過程とホワイトノイズ

ベクトル過程 $\boldsymbol{\varepsilon}_t$ が全ての時点 $t$ において次の二つの式を満たすとします。
$$
\begin{align}
E(\boldsymbol{\varepsilon}_t) &= \mathbf{0}&\\
E(\boldsymbol{\varepsilon}_t \boldsymbol{\varepsilon}_{t-k}^{\top}) &=
\left\{
\begin{matrix}
\boldsymbol{\Sigma},& k=0\\\
\mathbf{0},& k\neq 0 \
\end{matrix}
\right.
\end{align}
$$

この時、 $\boldsymbol{\varepsilon}_t$ をベクトルホワイトノイズと呼びます。
このベクトルホワイトノイズは弱定常であり、自己相関を持たないことはすぐわかります。
$k\neq 0$の時、$k$次の自己共分散行列が$\mathbf{0}$ですから。

沖本先生の本で注意されている通り、同時点での各変数は相関を持つことができます。

ただ、一点、次の記載があるのには少し引っかかります。(初版第13刷 P.76)

$\boldsymbol{\Sigma}$ は$n\times n$正定値行列であり、対角行列である必要はないことには注意が必要である.

$n\times n$正定値行列であることは確かなのでしが、これ、必ず対角行列になるような気がします。
というのも、$\boldsymbol{\Gamma}_k = \boldsymbol{\Gamma}^{\top}_{-k}$に、$k=0$ を代入したら、
$\boldsymbol{\Gamma}_0 = \boldsymbol{\Gamma}^{\top}_{0}$ となり、これは対角行列ですよね。

僕が何か勘違いしているのかもしれませんがおそらく誤植か何かのように思います。
(出版社サイトの今日時点の正誤表には入ってないようです)

弱定常ベクトル過程

ベクトル自己回帰モデルについて理解していくにあたって、弱定常過程やホワイトノイズの概念をベクトルに拡張したものを考える必要があるので、
ここで各用語の定義を整理しておきます。(前回の記事と順番前後してしまいました。)
出典は今回も沖本先生の本です。

また、1変数の場合については次の記事をご参照ください。
参考: 確率過程の弱定常性と強定常性

まず、 $\mathbf{y}_t = (y_{1t}, y_{2t}, \ldots, y_{nt})^{\top}$ は $n\times1$の行列(列ベクトル)としておきます。
そして、これの期待値ベクトルは普通に要素ごとの期待値をとって、次のように定義されます。
$$
E(\mathbf{y}_t) = (E(y_{1t}), \ldots, E(y_{nt}))^{\top}.
$$

次に、自己共分散を拡張します。ここからがポイントで、$y_{it}$に対して$y_{i,t-k}$との共分散を考えるだけではなく、
$\mathbf{y}_t$の他の要素の過去の値$y_{j,t-k}$との共分散も考慮します。そのため、結果は次のように行列になります。

そのため、$k$次自己共分散行列は次の$n\times n$行列で定義されます。

$$
\begin{align}
Cov(\mathbf{y}_t, \mathbf{y}_{t-k}) &= [Cov(y_{it},y_{j,t-k})]_{ij}\\
&=\left(\
\begin{matrix}\
Cov(y_{1t},y_{1,t-k}) & \cdots & Cov(y_{1t},y_{n,t-k}) \\\
\vdots & \ddots & \vdots \\\
Cov(y_{nt},y_{1,t-k}) & \cdots & Cov(y_{nt},y_{n,t-k}) \
\end{matrix}\
\right)
\end{align}.
$$

この行列の対角成分は、各変数の$k$次自己共分散に等しくなります。
また、$k=0$の時、$k$次自己共分散行列は、分散共分散行列に等しくなり、対称行列になりますが、
$k\neq 0$の時は、一般には対称行列にならないので注意が必要です。

期待値も$k$次自己共分散行列も定義式の中に$t$が入っているので、
これらは$t$の関数なのですが、この値がどちらも$t$に依存しなかった時、
1変数の場合と同じように、ベクトル過程は弱定常と言われます。

弱定常性を仮定する時は、期待値を$\boldsymbol{\mu}$、$k$次自己共分散行列を$\boldsymbol{\Gamma}_k$と書きます。
なお、$\boldsymbol{\Gamma}_k = \boldsymbol{\Gamma}^{\top}_{-k}$が成り立ちます。(転地を等号が成り立たない点に注意が必要です。)

さて、自己共分散行列は、異なる時点での各変数の関係を表したものですが、各要素の値の大小は単位に依存してしまい、
これだけ見ても関係が強いのか弱いのかわかりません。
そこで、自己共分散行列を標準化した自己相関行列を考えます。
それは次の式で定義されます。
$$
\boldsymbol{\rho}_k = Corr(\mathbf{y}_t, \mathbf{y}_{t-k}) = [Corr(y_{it},y_{j,t-k})]_{ij}.
$$

また、$\mathbf{D} = diag(Var(y_1), \ldots, Var(y_n))$、($\mathbf{y}_t$の分散を対角成分に持つ対角行列) とおくと、
次のように書けます。
$$
\boldsymbol{\rho}_k = \mathbf{D}^{-1/2} \boldsymbol{\Gamma}_k \mathbf{D}^{-1/2}.
$$

自己相関行列も自己共分散行列と同じように、対角成分は各変数の自己相関になり、
$\boldsymbol{\rho}_k = \boldsymbol{\rho}^{\top}_{-k}$ が成立します。

長くなってきたので一旦記事を切ります。(ベクトルホワイトノイズは次回。)

ベクトル自己回帰モデル

久しぶりに時系列分析の話題です。
今回はベクトル自己回帰(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歳を超えて未経験の状態からいきなりデータサイエンティストに転職してきた頃の話とか、
その後、分析チームのマネジメントや採用業務等も担当するようになっているので、
何かしら誰かの参考になる話もあるかと思います。