BLSTMを用いた文書分類でデザイナーズマンション予測に再挑戦

はじめに

仕事で深層学習を使うことはないのですが、扱える技術の幅を広げていきたいなと思い、BLSTMを用いた文書分類についてkerasでの簡単なサンプルコードをやってみようと思います。データは以前集めた、某不動産紹介サイトの賃貸マンションの設備に関する文書とそれがデザイナーズマンションかどうかのラベルを使います。そして文書内の単語からその文書がデザイナーズマンションかどうかを予測します。前回はAUCで83%だったので、それを超えれると良いですね。

目次

単純なRNNとは

  • モチベーション
    • フィードフォワード型のニューラルネットワークではうまく扱うことができない時系列データをうまく扱えるようにすること。
  •  特徴
    • 入力が互いに関係している(多層パーセプトロンの際は置かれていない仮定)
      • 直訳すると循環するニューラルネットワークとなる。
    • 最初の文の単語が2つ目、3つ目の単語に影響を与える可能性を考慮。
    具体的な関数の形としては、 $$ h_t = \tanh (W h_{t-1} + U x_t) \\\ y_t = softmax(V h_t) $$ で与えられる。
    \( h_t \)は隠れ層の状態を、\( x_t \)は入力変数を、\( y_t \)は出力ベクトルを表している。(他のアルファベットは重み行列)

LSTMとは

  • モチベーション
    • 単純なRNNの勾配消失問題を解決するために提案された手法。
  • 特徴
    • 単純なRNNに置き換えることで性能が大幅に向上することも珍しくない。
      • 時系列データ、長い文章、録音データからなる長期的なパターンを取り出すことを得意としている手法。
    • 勾配消失問題に強い
    • Long Short-Term Memory:長短期記憶
      • 長期依存性を学習できるRNNの亜種。
      • 入力ゲート、忘却ゲート、出力ゲート、内部隠れ層という4つの層が相互に関わり合う。重要な入力を認識し、それを長期状態に格納すること、必要な限りで記憶を保持すること、必要なときに記憶を取り出すことを学習する。
        • 入力ゲート:内部隠れ層のどの部分を長期状態に加えるかを決める。
        • 忘却ゲート:長期状態のどの部分を消去するか決める。
        • 出力ゲート:各タイムステップで、長期状態のどの部分を読み出し、出力するかを決める。
$$
i = \sigma (W_i h_{t-1} + U_i x_t) \\
f = \sigma (W_f h_{t-1} + U_f x_t) \\
o = \sigma (W_o h_{t-1} + U_ox_t ) \\
g = \tanh (W_g h_{t-1} + U_g x_t) \\
c_t = (c_{t-1} \otimes f ) \otimes ( g \otimes i) \\
h_t = \tanh (c_t) \otimes o
$$ i:入力ゲート(input)
f:忘却ゲート(forget)
o:出力ゲート(output)
\( \sigma \):シグモイド関数
g:内部隠れ層状態
\(c_t \):時刻tにおけるセル状態
\(h_t \):時刻tにおける隠れ状態

Bidirectional LSTMとは

  • モチベーション
    •  従来のRNNの制約を緩和するために導入。
  • 特徴
    • ある特定の時点で過去と将来の利用可能な入力情報を用いて訓練するネットワーク
      • 従来のRNNのニューロンをフォワード(未来)なものとバックワード(過去)なものとの2つに分ける。
        • 2つのLSTMを訓練している。
      • 全体のコンテキストを利用して推定することができるのが良いらしい。
        • 文章で言うと、文章の前(過去)と後(未来)を考慮して分類などを行うイメージ。
      • 人間もコンテキストを先読みしながら解釈することもある。
    • BLSTMは全てのシークエンスの予測問題において有効ではないが、適切に扱えるドメインで良い結果を残す。
    • 1997年とかなり歴史があるものらしい
出典:Deep Dive into Bidirectional LSTM

Bidirectional LSTMで文書分類

今回は、BLSTMを使って、マンションの設備に関するテキスト情報から、そのマンションがデザイナーズマンションかどうかを予測します。全体のソースコードはGoogle Colabを御覧ください。

データ

データを確認してみると、
 
今回のデータがテキストとラベルからなっていることがわかる。なお、データ数は1864件あり、そのうち16%がデザイナーズマンションです。

前処理

Google Colab上で形態素解析を行うために、MeCabをインストールする。
 
これでMaCab NeologdをGoogle ColabのPythonで実行できます。
テキストデータの前処理を行うために名詞だけを抽出してリスト形式で返す関数を定義します。
 
ここで、語彙に関するデータを作成します。
 
続いて、単語とそのインデックスからなるdict形式のデータを作成します。
 
最後に、これまでに作成した語彙とそのインデックスのデータを用いて、実際に学習に使う入力データと出力データを作成します。
 

実行

先程作成したデータを訓練データとテストデータに分けます。
ここで、AUCを計算するための関数を定義します。(まるまる拝借しました)
続いて、kerasを用いて、ネットワークを構築し、コンパイルを行います。
  テストデータでの予測精度を確認します。
実際のラベルとテキストデータを見比べてみます。
LSTMの場合は、BLSTMに関する記述(1行)をネットワークから除外すれば良いので簡単に試すことができます。 せっかくなので比較してみたところ、AUCは0.841でした。BLSTMは0.844だったので、若干上回る程度のようです。 以前扱った記事ではデザイナーズマンション分類はAUCが0.83程度だったので、LSTMを使ったほうが精度が少しだけ高いようです。

追記

TensorBoardの検証のloglossを見てみます。 エポックに対してloglossが非常に荒れています。学習曲線としてはセオリー的にアウトなようです。一方で、検証のAUCはエポックに従って高まるようです。 AUCは良くなっているけどloglossが増え続けている。loglossが下がらないと過学習している可能性が高いので、これは過学習しているだけなのだろう。

参考文献

 

[数理統計学]正規分布から導かれる分布(カイ二乗分布/t分布/F分布)の期待値と分散の導出まとめ

はじめに

前回は、連続型確率分布に関するよくある分布とその平均と分散の導出についてひたすら記しました。今回は同じく連続型ながら、正規分布から導かれる確率分布に関して同様に記していこうと思います。大学のときに数理統計学の授業でがっつりやったはずの領域だったのですが、見事に忘れていたので、電車の中などでしっかりと復習していこうと思います。

※PC版でないと数式が切れてしまうので、SP版での閲覧はおすすめしません。私、SKUEのスマホだけは表示崩れがないようにチェックしております。通勤電車で定期的に読むので。

今回登場する連続分布

『数理統計学―基礎から学ぶデータ解析』という本に登場するものを式を補いながら紹介します。

1.カイ二乗分布

ガンマ分布\( \Gamma \left ( \frac{n}{2}, \frac{1}{2} \right ) \)、すなわち、密度関数が以下の式で与えられるとき、この確率分布を自由度nの\( \chi^2 \)分布という。

$$ f(x ; n) = \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}-1} e^{-\frac{x}{2}} \\
x > 0 \\
n = 1, 2, \dots
$$

\( \chi^2 \)分布の再生性

XとYがカイ二乗分布に従うとして、それらの和がカイ二乗分布になる性質を再生性と呼ぶ。
$$ X \sim \chi^2 (m) \\
Y \sim \chi^2 (n) \\
X + Y \sim \chi^2 (m+n)
$$

以下ではこの再生性に関して証明を行う。
ここで\( z = X + Y \)とする。zの密度関数は、以下で表される。

$$ g(z) = \int_0^z f(x ; m) f(y ; n) dx \\
= \int_0^z f(x ; m) f(z-x ; n) dx \\
= \int_0^z \frac{ x^{\frac{m}{2}-1} e^{-\frac{x}{2}} }{2^{\frac{m}{2}} \Gamma \left ( \frac{m}{2} \right ) } \frac{ (z-x)^{\frac{n}{2}-1} e^{-\frac{z-x}{2}} }{2^{\frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) } dx
$$

ここで、\( u = \frac{x}{z} \)として、変数変換を行う。( \( \frac{dx}{du} = z \quad , x = zu \) )

$$ = \int_0^1 \frac{ (zu)^{\frac{m}{2}-1} e^{-\frac{zu}{2}} }{2^{\frac{m}{2}} \Gamma \left ( \frac{m}{2} \right ) } \frac{ (z-zu)^{\frac{n}{2}-1} e^{-\frac{z(1-u)}{2}} }{2^{\frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) } z du $$

$$ = \frac{z^{\frac{m}{2} + \frac{n}{2} -1} e^{-\frac{z}{2}} }{ 2^{\frac{m}{2} + \frac{n}{2} } \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \int_0^1 u^{\frac{u}{2}-1}(1-u)^{\frac{u}{2}-1}du $$

積分記号以降の部分はベータ分布であるから、

$$ = \frac{ e^{-\frac{z}{2}} z^{\frac{m+n}{2} -1} }{ 2^{\frac{m+n}{2} } \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } B \left ( \frac{m}{2}, \frac{n}{2} \right ) $$

以前記したベータ分布とガンマ関数との関係( \( B(x, y) = \frac{\Gamma(x) \Gamma(y)}{\Gamma(x+y)} \) )より、

$$ = \frac{ e^{-\frac{z}{2}} z^{\frac{m+n}{2} -1} }{ \Gamma \left ( \frac{m+n}{2} \right ) 2^{\frac{m+n}{2} } } $$

これは\( \chi^2 (m+n) \)の確率密度関数である。

定理:\( X_1, X_2, \dots , X_n \)は独立に正規分布\( N(0, 1) \)に従うとき、\( Y = X_1^2 + X_2^2 + \dots + X_n^2 \)は\( \chi^2(n) \)に従う。

以下はUCLAのレクチャーノートに従っています。

まず、\( X_i^2 \)がカイ二乗分布に従うことを示します。
標準正規分布に従う変数zの二乗としてxを定義します。
$$ x = z^2 \\
f(z) = \frac{1}{\sqrt {2 \pi } } e^{-\frac{1}{2}z^2}
$$
xの確率は以下のように表される。
$$ F_X (x) = P( X \le x ) \\
= P(Z^2 \le x ) \\
= P( – \sqrt {x} \le Z \le \sqrt {x} )
$$

この場合、求めたい確率は1から両端を差し引いたものとなる。

$$ {left} = F_z(-\sqrt {x}) \\
{right} = 1 – F_z(\sqrt {x}) \\
F_X (x) = 1 – {left} – {right} \\
= F_z(\sqrt {x}) – F_z(-\sqrt {x})
$$

確率密度を求めるために両辺をxで微分し整理する。

$$ \frac{d F_X (x) }{dx} = \frac{d F_z(\sqrt {x}) }{ d \sqrt {x} } \frac{d \sqrt {x}}{dx} – \frac{d F_z( – \sqrt {x}) }{ d – \sqrt {x} } \frac{d – \sqrt {x}}{dx} \\\\
= \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}x}\frac{1}{2}x^{-\frac{1}{2}} – (-1) \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}x}\frac{1}{2}x^{-\frac{1}{2}} \\\\
= \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}x}x^{-\frac{1}{2}} \\\\
= \frac{x^{-\frac{1}{2}} e^{-\frac{x}{2}}}{2^{\frac{1}{2}} \Gamma \left ( \frac{1}{2} \right ) } \sim \chi_1^2
$$

これは自由度1のカイ二乗分布であり、正規分布に従う変数の二乗は自由度1のカイ二乗分布に従うことが示された。

ここで、ガンマ分布の積率母関数について導出を行う。

$$ M_{X}(t) = E \left ( e^{tX} \right ) \\
= \int _{0}^{\infty} e^{tx} \frac{x^{\alpha – 1} e^{-\frac{x}{\beta}}}{\beta^{\alpha} \Gamma \left ( \alpha \right )} dx \\
= \frac{1}{\beta^{\alpha} \Gamma \left ( \alpha \right )}\int _{0}^{\infty} x^{\alpha – 1} e^{-x \left ( \frac{1-\beta t}{\beta} \right ) }dx
$$

ここで以下の変数変換を行う。
$$ y = x \left ( \frac{1-\beta t}{\beta} \right ) $$
つまり、\( x = \frac{\beta}{1-\beta t} y \)、\( \frac{dx}{dy} = \frac{\beta}{1-\beta t} \)であるから、これらを代入すると、

$$ M_X(t) = \frac{1}{\beta^{\alpha} \Gamma \left ( \alpha \right )} \int _{0}^{\infty} \left ( \frac{\beta}{1-\beta t} \right )^{\alpha – 1}y^{\alpha – 1}e^{-y}\frac{\beta}{1-\beta t}dy \\\\
= \frac{1}{\beta^{\alpha} \Gamma \left ( \alpha \right )} \left ( \frac{\beta}{1-\beta t} \right )^{\alpha – 1} \frac{\beta}{1-\beta t} \int _{0}^{\infty} y^{\alpha – 1}e^{-y} dy
\\\\
= (1-\beta t)^{-\alpha}
$$

よって、先程求めた自由度1のカイ二乗分布は積率母関数で表すと以下のようになる。
$$ M_X(t) = M_{Z^2}(t) = (1-2t)^{-\frac{1}{2}} $$

今回の定理において、\( X_i^2 \)は各々が独立なので、先程求めた自由度1のカイ二乗分布の積率母関数を掛け合わせることで\( Y = \sum_{i}^{n} X_i^2 \)の積率母関数は求まる。

$$ M_Y(t) = M_{Z_1^2}(t) \times M_{Z_2^2}(t) \times \dots \times M_{Z_n^2}(t) \\\\
= (1-2t)^{-\frac{1}{2}} \times (1-2t)^{-\frac{1}{2}} \times \dots \times (1-2t)^{-\frac{1}{2}} \\\\
= (1-2t)^{-\frac{n}{2}}
$$

これは自由度nのカイ二乗分布の積率母関数と同じであるから、\( Y = \sum_{i}^{n} X_i^2 \)が自由度nのカイ二乗分布に従うことが示された。

カイ二乗分布の平均

$$ E \left ( X \right ) = \int _{0}^{\infty} x \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx \\
= \int _{0}^{\infty} \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}} e^{-\frac{x}{2}}dx
$$

ここで、部分積分の公式より

$$ \int _{0}^{\infty} x^{\frac{n}{2}} e^{-\frac{x}{2}}dx \\\\
= \left [ \frac{2}{n+2}x^{\frac{n}{2}+1} e^{-\frac{x}{2}} \right ]_0^{\infty}- \frac{n}{2}(-2) \int _{0}^{\infty} x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx
\\\\
= n \int _{0}^{\infty} x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx
$$

となるから、カイ二乗分布の積分が1になることから、平均は

$$ E \left ( X \right ) = n \int _{0}^{\infty} \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx \\
= n
$$
となる。

カイ二乗分布の分散

分散は
$$ V(X) = E \left ( X^2 \right ) – \left [ E(X) \right ]^2 $$
より求まるので、

$$ E(X^2) = \int _{0}^{\infty} x^2 \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx \\
= \int _{0}^{\infty} \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}+1} e^{-\frac{x}{2}}dx
$$

ここで、部分積分の公式より

$$ \int _{0}^{\infty} x^{\frac{n}{2}+1} e^{-\frac{x}{2}}dx \\\\
= \left [ \frac{2}{n+4}x^{\frac{n}{2}+2} e^{-\frac{x}{2}} \right ]_0^{\infty}- \frac{n+2}{2}(-2) \int _{0}^{\infty} x^{\frac{n}{2}} e^{-\frac{x}{2}}dx \\\\
= (n+2)\int _{0}^{\infty} x^{\frac{n}{2}} e^{-\frac{x}{2}}dx \\\\
= (n+2) \left [ \left [ \frac{2}{n+2}x^{\frac{n}{2}+1} e^{-\frac{x}{2}} \right ]_0^{\infty}- \frac{n}{2}(-2) \int _{0}^{\infty} x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx \right ] \\\\
= n(n+2) \int _{0}^{\infty} x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx
$$

となるから、カイ二乗分布の積分が1になることから、

$$ E(X^2) = n(n+2) \int _{0}^{\infty} \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx \\\\
= n(n+2)
$$

よってカイ二乗分布の分散は
$$ V(X) = n(n+2) – n^2 \\
= 2n $$
となる。

2.t分布

以下の確率密度関数で与えられる分布をt分布と呼ぶ。

$$ f(x ; n) = \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} } \\
-\infty < x < \infty \\
n = 1, 2, \dots
$$

確率密度関数から明らかなように、パラメータは自由度nだけからなる。
なお、自由度が1のとき、t分布はコーシー分布と同じになる。
$$ f(x ; 1) = \frac{1}{\pi} \frac{1}{1 + x^2} $$

t分布の特徴としては、
・y軸に関して対称
・標準正規分布よりも左右に裾が重い
・自由度が1よりも大きいとき、xf(x)は可積分となり、平均は0となる。
・自由度が1の場合はコーシー分布と同じであるため平均も分散も存在しない。
などがあげられる。

定理:Xは\( N(0, 1) \)に、Yは\( \chi^2(n) \)にそれぞれ従い、かつ互いに独立とする。このとき、\( t = \frac{x}{\sqrt {y/n}} \)は\( t(n) \)に従う。

XとYの結合密度関数は
$$ f(x, y) = \frac{1}{\sqrt {2\pi}} e^{ -\frac{x^2}{2}} \frac{ y^{ \frac{n}{2} -1} }{ 2^{ \frac{n}{2} } \Gamma \left ( \frac{n}{2} \right ) } e^{ – \frac{y}{2} } $$
となる。
$$ t = \frac{x}{\sqrt {y / n} } \\
x = t \sqrt {y / n} = t \sqrt {s / n} \\
y = s
$$
とおいて、確率の変数変換を行うと、tとsの結合密度関数は、

$$ g(t, s) = f \left ( t \sqrt {s / n}, s \right ) | J | \\
= f \left ( t \sqrt {s / n}, s \right ) \begin{bmatrix} \frac{\partial x}{\partial t} & \frac{\partial x}{\partial s} \\ \frac{\partial y}{\partial t} & \frac{\partial y}{\partial s} \end{bmatrix} \\
= f \left ( t \sqrt {s / n}, s \right ) \sqrt {s / n} \\
= \frac{1}{\sqrt {2\pi}} e^{ -\frac{t^2 \frac{s}{n} }{2}} \frac{ s^{ \frac{n}{2} -1} }{ 2^{ \frac{n}{2} } \Gamma \left ( \frac{n}{2} \right ) } e^{ – \frac{s}{2} } \sqrt {s / n} \\
= \frac{ s^{ \frac{1}{2}(n-1)} e^{-\frac{1}{2} \left ( 1 + \frac{t^2}{n}s \right ) } }{ \sqrt {2\pi} 2^{ \frac{n}{2} } \Gamma \left ( \frac{n}{2} \right ) \sqrt {n} }
$$

ここで\( u = \frac{1}{2} \left ( 1 + \frac{t^2}{n} \right ) s \)とおく。なお、\( s = \frac{2u}{ \left ( 1 + \frac{t^2}{n} \right ) } \)、\( \frac{ds}{du} = \frac{2}{1 + \frac{t^2}{n}} \)であるからこれらを代入し、分子について定積分をとりtを固定し、sを変換する。

$$ \int _{0}^{\infty} s^{ \frac{1}{2} (n-1)} e^{ – \frac{1}{2} \left ( 1 + \frac{t^2}{n} \right )s }ds \\
= \int _{0}^{\infty} \left ( \frac{2u}{ 1 + \frac{t^2}{n} } \right )^{ \frac{1}{2}(n-1) }e^{-u} \frac{ds}{du}du \\
= \int _{0}^{\infty} \left ( \frac{2u}{ 1 + \frac{t^2}{n} } \right )^{ \frac{1}{2}(n-1) }e^{-u} \frac{2}{1 + \frac{t^2}{n}} du \\
= \frac{ 2^{ \frac{1}{2} (n-1) + 1 } }{ \left ( 1 + \frac{t^2}{n} \right )^{ \frac{1}{2}(n-1) + 1 } } \int _{0}^{\infty} u^{ \frac{1}{2}(n-1)} e^{-u}du \\
= \frac{ 2^{ \frac{1}{2} (n+1) } }{ \left ( 1 + \frac{t^2}{n} \right )^{ \frac{1}{2}(n+1)} } \int _{0}^{\infty} u^{ \frac{1}{2}(n-1)} e^{-u}du
$$

ここで、ガンマ関数の定義より、
$$ \Gamma ( \alpha ) = \int_{0}^{\infty} u^{ \alpha – 1} e^{-u}du \\
\alpha – 1 = \frac{1}{2}(n-1) \\
\alpha = \frac{1}{2}(n+1)
$$
であるから、これを考慮すると、
$$ = \frac{ 2^{ \frac{1}{2} (n+1) } }{ \left ( 1 + \frac{t^2}{n} \right )^{ \frac{1}{2}(n+1)} } \Gamma \left ( \frac{1}{2} (n+1) \right ) $$
となる。

以上より、tの密度関数は

$$ g(t) = \int_{0}^{\infty} g(t, s) ds \\\\
= \frac{ 2^{ \frac{1}{2} (n+1) } }{ \left ( 1 + \frac{t^2}{n} \right )^{ \frac{1}{2}(n+1)} } \Gamma \left ( \frac{1}{2} (n+1) \right ) \frac{1}{\sqrt {2\pi} 2^{ \frac{n}{2} } \Gamma \left ( \frac{n}{2} \right ) \sqrt {n}} \\\\
= \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt{n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{t^2}{n} \right )^{ \frac{n+1}{2}} }
$$

これはt分布の密度関数と同じである。

t分布の平均

$$ E \left ( X \right ) = \int _{-\infty}^{\infty} x \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} } dx $$
奇関数であることから
$$ = 0 $$
となる。

t分布の分散

分散は
$$ V(X) = E \left ( X^2 \right ) – \left [ E(X) \right ]^2 \\
= E \left ( X^2 \right ) $$
より求まる。

$$ E \left ( X^2 \right ) = \int _{-\infty}^{\infty} x^2 \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} } dx \\\\
= \int _{-\infty}^{0} x^2 \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} } dx + \int _{0}^{\infty} x^2 \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} } dx
$$

偶関数であることから以下のように表される。
$$ = 2 \int _{0}^{\infty} x^2 \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} } dx \\
= 2 \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right )} \int _{0}^{\infty} \frac{x^2}{ \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} }dx
$$

ここで、\( y = \frac{1}{1+\frac{x^2}{n}} \)により変数変換を行う。
なお、\( \left ( 1 + \frac{x^2}{n} \right ) = \frac{1}{y} \)、\( x = \left [ \frac{(1-y)n}{y} \right ] ^{\frac{1}{2}} \)、
\( \frac{dx}{dy} = \frac{ (1-y)^{-\frac{1}{2}} n^{\frac{1}{2}}(-1)y^{-\frac{3}{2}}}{2} \)となる。

$$ = (-1)2\frac{ \Gamma \left ( \frac{n+1}{2} \right ) n^{\frac{3}{2}} }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right )} \int _{0}^{\infty} \frac{(1-y)^{\frac{1}{2}} y^{\frac{n+1}{2}-1-\frac{3}{2}} }{2}dy \\\\
= 2\frac{ \Gamma \left ( \frac{n+1}{2} \right ) n^{\frac{3}{2}} }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right )} \int _{0}^{\infty} \frac{(1-y)^{\frac{1}{2}} y^{\frac{n}{2}-2} }{2}dy
$$

ベータ関数( \( B( \alpha, \beta ) = \int_0^{1} t^{\alpha – 1}(1-t)^{ \beta – 1} dt \) )より、

$$ = \frac{ \Gamma \left ( \frac{n+1}{2} \right ) n^{\frac{3}{2}} }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right )} B\left ( \frac{n-2}{2}, \frac{3}{2} \right ) $$

ベータ関数とガンマ関数の関係より、
$$ B\left ( \frac{n-2}{2}, \frac{3}{2} \right ) = \frac{ \Gamma \left ( \frac{n-2}{2} \right ) \Gamma \left ( \frac{3}{2} \right ) }{ \Gamma \left ( \frac{n-2}{2} + \frac{3}{2} \right ) } \frac{1}{2} \\
= \frac{ \Gamma \left ( \frac{n-2}{2} \right ) \Gamma \left ( \frac{3}{2} \right ) }{ \Gamma \left ( \frac{n}{2} + \frac{1}{2} \right ) }
$$

よって、
$$ = \frac{ \Gamma \left ( \frac{n+1}{2} \right ) n^{\frac{3}{2}} }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right )} \frac{ \Gamma \left ( \frac{n-2}{2} \right ) \Gamma \left ( \frac{3}{2} \right ) }{ \Gamma \left ( \frac{n}{2} + \frac{1}{2} \right ) } $$

ガンマ関数は\( \Gamma (n) = (n-1) \Gamma (n-1) \)となることから、\( \Gamma (n+1) = n \Gamma (n) \)あるいは\( \Gamma \left( \frac{n}{2} – 1 \right) = \frac{ \Gamma \left ( \frac{n}{2}
\right ) }{\frac{n}{2} – 1} \)であることを用いると、
$$ = \frac{n^{ \frac{3}{2}} \Gamma \left ( \frac{n}{2} -1 \right ) \Gamma \left ( \frac{n}{2} + 1 \right ) }{ n^{\frac{1}{2}} \sqrt{\pi} \Gamma \left ( \frac{n}{2} \right ) } \\
= \frac{n \Gamma \left ( \frac{n}{2} \right ) }{ \pi^{\frac{1}{2}} \Gamma \left ( \frac{n}{2} \right ) \left ( \frac{n}{2} – 1 \right )} \frac{1}{2} \Gamma \left ( \frac{n}{2} \right ) \\
= \frac{n}{n-2}
$$
となる。

よって、
$$ V(X) = \frac{n}{n-2} $$

3.F分布

以下の確率密度関数で与えられる分布をF分布と呼ぶ。
$$ f(x;m;n) = \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ x^{ \frac{m}{2} -1}}{ (mx+n)^{ \frac{m+n}{2} } } \\
x > 0 \\
m, n = 1, 2, \dots
$$

定理:Xは\( \chi^2(m) \)に、Yは\( \chi^2(n) \)に独立に従うとする。このとき、\( F = \frac{ \frac{1}{m}X }{ \frac{1}{n}Y } \)は\( F(m, n) \)に従う。

ここで、\( (X, Y) \rightarrow (F, U) \)と確率の変数変換を行う。

\(
\begin{cases}
F = \frac{ \frac{1}{m}X }{ \frac{1}{n}Y } \\
U = Y \\
\end{cases}
\)
また、\( X = \frac{1}{n}YFM \)、\( Y = U \)となる。

よって確率の変数変換の際のヤコビアンは、

$$ \begin{bmatrix} \frac{\partial X}{\partial F} & \frac{\partial X}{\partial U} \\ \frac{\partial Y}{\partial F} & \frac{\partial Y}{\partial U} \end{bmatrix} = \begin{bmatrix} \frac{1}{n}Um & 0 \\ 0 & 1 \end{bmatrix} = \frac{m}{n}U $$
となる。

FとUの結合分布の密度関数は
$$ h(x, y) = \frac{ x^{ \frac{m}{2} -1} e^{ – \frac{x}{2} }} { 2^{ \frac{m}{2}} \Gamma \left ( \frac{m}{2} \right )} \frac{ y^{ \frac{n}{2} -1} e^{ – \frac{y}{2} }} { 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right )} $$
となる。
確率の変数変換を行い整理すると、
$$ g(f,u) = h(x, y) |J | \\
= h \left ( \frac{1}{n}ufm , u \right ) \left | \frac{m}{n}u \right | \\
= \frac{ u^{\frac{m+n}{2} -1} n^{-\frac{m}{2}} m^{\frac{m}{2}} f^{\frac{m}{2} -1} e^{-\frac{u}{2} \left ( \frac{mf}{n} + 1 \right ) } }{ 2^{ \frac{m}{2} } 2^{ \frac{n}{2} } \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) }
$$

ここでfを固定して、\( t = \frac{1}{2} \left ( \frac{mf}{n} + 1 \right ) u \)として変数変換を行う。なお、\( u = \frac{2t}{ \frac{mf}{n}+1 } \)、\( \frac{du}{dt} = \frac{2}{ \frac{mf}{n} + 1 } \)となる。

$$ g(f) = \int_0^{\infty} g(f, u)du \\\\
= \frac{ n^{-\frac{m}{2}} m^{\frac{m}{2}} f^{\frac{m}{2} -1} }{ 2^{ \frac{m}{2} } 2^{ \frac{n}{2} } \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \int_0^{\infty} e^{-\frac{u}{2} \left ( \frac{mf}{n} + 1 \right ) } u^{\frac{m+n}{2} -1} du
$$

ここで積分の中についてのみ注目すると、

$$ \int_0^{\infty} e^{-\frac{u}{2} \left ( \frac{mf}{n} + 1 \right ) } u^{\frac{m+n}{2} -1} du \\\\
= \int_0^{\infty} e^{-t} \frac{2}{ \frac{mf}{n} + 1 } \left ( \frac{2t}{ \frac{mf}{n} + 1} \right )^{ \frac{m+n}{2}-1 }dt \\\\
= \left ( \frac{2n}{mf+n} \right )^{\frac{m+n}{2}} \int_0^{\infty} t^{ \frac{m+n}{2}-1}e^{-t}dt
$$

よって、

$$ g(f) = \frac{ n^{-\frac{m}{2}} m^{\frac{m}{2}} f^{\frac{m}{2} -1} }{ 2^{ \frac{m}{2} } 2^{ \frac{n}{2} } \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \left ( \frac{2n}{mf+n} \right )^{\frac{m+n}{2}} \int_0^{\infty} t^{ \frac{m+n}{2}-1}e^{-t}dt $$

ガンマ関数の定義より、
$$ \Gamma ( \alpha ) = \int_{0}^{\infty} u^{ \alpha – 1} e^{-u}du \\
\alpha – 1 = \frac{m+n}{2} – 1 \\
\alpha = \frac{m+n}{2}
$$
であるから、これを考慮し整理する。

$$ g(f) = \frac{ n^{-\frac{m}{2}} m^{\frac{m}{2}} f^{\frac{m}{2} -1} }{ 2^{ \frac{m}{2} } 2^{ \frac{n}{2} } \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \left ( \frac{2n}{mf+n} \right )^{\frac{m+n}{2}} \Gamma \left ( \frac{m+n}{2} \right ) \\\\
= \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ f^{ \frac{m}{2} -1}}{ (mf+n)^{ \frac{m+n}{2} } }
$$

これはF分布の確率密度関数である。

F分布の平均

$$ E(X) = \int_0^{\infty} x \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ x^{ \frac{m}{2} -1}}{ (mx+n)^{ \frac{m+n}{2} } } dx \\\\
= \int_0^{\infty} \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ x^{ \frac{m}{2}}}{ (mx+n)^{ \frac{m+n}{2} } } dx
$$

ここで、\( y = \frac{x}{mx+n} \)と変数変換すると、\( x = \frac{ny}{1-my} \)、\( \frac{dx}{dy} = \frac{n}{(1-my)^2} \)、\( mx + n = \frac{nym}{1-my} + n = \frac{n}{1-my} \)であるから、これらを代入すると、

$$ E(X) = \int_0^{\infty} \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ \left ( \frac{ny}{1-my} \right )^{\frac{m}{2}} }{ \left ( \frac{n}{1-my} \right )^{\frac{m+n}{2}} } \frac{n^2}{n (1-my)^2 }dy \\\\
= \int_0^{\infty} \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } y^{\frac{m+n}{2}} \left ( \frac{ny}{1-my} \right )^{-\frac{n}{2}} \frac{1}{n} \left ( \frac{n}{1-my} \right )^2 dy
$$

ここで、\( z=my \)と変数変換すると、\( y = \frac{z}{m} \)、\( \frac{dy}{dz} = \frac{1}{m} \)であるから、これらを代入すると、

$$ E(X) = \int_0^{\infty} \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \left ( \frac{z}{m} \right )^{\frac{m+n}{2}} \left ( \frac{ny}{1-z} \right )^{-\frac{n}{2}} \frac{1}{n} \left ( \frac{n}{1-z} \right )^{2}\frac{1}{m}dz \\\\
= \frac{ \Gamma \left ( \frac{m+n}{2} \right ) }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) }m^{-1}n \int_0^{\infty} z^{ \frac{m}{2}}(1-z)^{\frac{n}{2}-2}dz
$$

ベータ関数、\( B(\alpha , \beta) = \int_0^{1} t^{\alpha – 1}(1-t)^{\beta -1}dt \)の定義より、\( \alpha = \frac{m}{2} + 1 \)、\( \beta = \frac{n}{2} -1 \)であるから、
$$ B \left ( \frac{m}{2}+1, \frac{n}{2} -1 \right ) = \int_0^{\infty} z^{ \frac{m}{2}}(1-z)^{\frac{n}{2}-2}dz $$
ベータ関数とガンマ関数の関係式より、
$$ B \left ( \frac{m}{2}+1, \frac{n}{2} -1 \right ) = \frac{ \Gamma \left ( \frac{m}{2} + 1 \right ) \Gamma \left ( \frac{n}{2} – 1 \right ) }{ \Gamma \left ( \frac{m+n}{2} \right ) } $$
となる。
よって、
$$ E(X) \\
= \frac{ \Gamma \left ( \frac{m+n}{2} \right ) }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) }m^{-1}n \frac{ \Gamma \left ( \frac{m}{2} + 1 \right ) \Gamma \left ( \frac{n}{2} – 1 \right ) }{ \Gamma \left ( \frac{m+n}{2} \right ) } \\
= \frac{n}{n-2}
$$
となる。

F分布の分散

分散は
$$ V(X) = E \left ( X^2 \right ) – \left [ E(X) \right ]^2 $$
より求まる。

$$ E(X^2) = \int_0^{\infty} x^2 \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ x^{ \frac{m}{2} -1}}{ (mx+n)^{ \frac{m+n}{2} } } dx \\\\
= \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \int_0^{\infty} \frac{ x^{ \frac{m}{2} -1}}{ (mx+n)^{ \frac{m+n}{2} } } dx
$$

ここで、先ほどと同様に変数変換を行うが、ここでは積分以降にのみ注目する。

$$ \int_0^{\infty} \frac{ x^{ \frac{m}{2} -1}}{ (mx+n)^{ \frac{m+n}{2} } } dx \\\\
= \int_0^{\infty} \left ( \frac{ny}{1-my} \right )^{ \frac{m}{2}+1 } \left ( \frac{n}{1-my} \right )^{ -\frac{m+n}{2}}\frac{1}{n} \left ( \frac{n}{1-my} \right )^{2} dy \\\\
= \int_0^{\infty} \left ( \frac{ny}{1-my} \right )^{ -\frac{n}{2}+1 } \frac{1}{n} \left ( \frac{n}{1-my} \right )^{2} dy
$$

ここで、\( z = my \)と変数変換を行う。なお、\( y = \frac{z}{m} \)、\( \frac{dy}{dz} = \frac{1}{m} \)であるから、これらを代入すると、

$$ = \int_0^{\infty} \left ( \frac{z}{m} \right )^{ \frac{m+n}{2}} \left ( \frac{nz/m}{1-z} \right )^{ – \frac{n}{2} +1} \frac{1}{n} \left ( \frac{n}{1-z} \right )^{2} \frac{1}{m} dz \\\\
= m^{ -\frac{m}{2} -2 }n^{- \frac{n}{2} + 2} \int_0^{\infty} z^{ \frac{m}{2} +1 }(1-z)^{ \frac{n}{2} -3 }dz
$$

ベータ関数、\( B(\alpha , \beta) = \int_0^{1} t^{\alpha – 1}(1-t)^{\beta -1}dt \)の定義より、\( \alpha = \frac{m}{2} + 2 \)、\( \beta = \frac{n}{2} -2 \)であるから、

$$ B \left ( \frac{m}{2}+2, \frac{n}{2} -2 \right ) = \int_0^{\infty} z^{ \frac{m}{2} +1 }(1-z)^{ \frac{n}{2} -3 }dz $$

よって、

$$ E(X^2) = \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } m^{ -\frac{m}{2} -2 }n^{- \frac{n}{2} + 2} B \left ( \frac{m}{2}+2, \frac{n}{2} -2 \right ) $$

ベータ関数とガンマ関数の関係式より、
$$ B \left ( \frac{m}{2}+2, \frac{n}{2} -2 \right ) = \frac{ \Gamma \left ( \frac{m}{2} + 2 \right ) \Gamma \left ( \frac{n}{2} – 2 \right ) }{ \Gamma \left ( \frac{m+n}{2} \right ) } $$
となる。
よって、

$$ E(X^2) = \frac{n^2}{m^2} \frac{ \Gamma \left ( \frac{m+n}{2} \right ) }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ \Gamma \left ( \frac{m}{2} + 2 \right ) \Gamma \left ( \frac{n}{2} – 2 \right ) }{ \Gamma \left ( \frac{m+n}{2} \right ) } \\\\
= \frac{n^2}{m^2} \frac{ \Gamma \left ( \frac{m+n}{2} \right ) }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ \left ( \frac{m}{2} + 2 \right ) \frac{m}{2} \Gamma \left ( \frac{m}{2} \right ) }{ \Gamma \left ( \frac{m+n}{2} \right ) } \frac{ \Gamma \left ( \frac{n}{2} \right ) }{ \left ( \frac{n}{2} – 2 \right ) \left ( \frac{n}{2} -1 \right ) } \\\\
= \frac{n^2 \left ( \frac{m+2}{2} \right ) \frac{m}{2} }{m^2 \left ( \frac{n-4}{2} \right ) \left ( \frac{n-2}{2} \right ) } \\\\
= \frac{n^2(m+2)}{m(n-4)(n-2)}
$$

分散の式にあてはめると、
$$ V(X) = E \left ( X^2 \right ) – \left [ E(X) \right ]^2 \\
= \frac{n^2(m+2)}{m(n-4)(n-2)} – \frac{n^2}{(n-2)^2} \\
= \frac{2n^2(n+m-2)}{m(n-4)(n-2)^2}
$$
となる。

参考文献

[1] Distributions related to the normal distribution
[2] 鈴木・山田 (1996), 『数理統計学―基礎から学ぶデータ解析』, 内田老鶴圃

[R]ボージョレ・ヌーボーのコメントに対してLDATSパッケージを使って時系列トピックモデルを扱う

はじめに

先日、某勉強会でLTをしました。その際に10秒だけ紹介したRのパッケージについて記事を書いてみようと思います。

LDATSパッケージについて

時系列でのトピックモデルを推定することができるパッケージです。
やっていることとしてはLDAでトピックを推定して次元を減らし、そのトピックの多変量時系列に関してベイズ手法による変化点検知のためのパラメータ推定を行っているようです。GitHubの該当しそうなソースコードに多変量のデータに対するsoftmax関数での回帰をやっているとの記述がある。(multinomial Bayesian Time Series analysis

元となっている論文を見る限り、BoW(Bag of Words)を想定して作っておらず、20~30程度のグループからなるデータに対して適用するのがちょうど良いです。アクセスログのページカテゴリや、マーケティングの顧客セグメントであればそんなに数は多くないので扱いやすいと思います。

データ

Webサイトから集めてきたボージョレ・ヌーボーのキャッチコピー14年分を今回は扱います。実は販売店側のキャッチコピーとワイン委員会が決めた評価が存在します。私の知っている世界は販売店側のキャッチコピーだけでした。

試してみた

今回はとにかく動くことだけを考えて、汚いコードとなっております。やっていることとしては、キャッチコピーを販売側とワイン委員会側のものを一つにつないで、数字を正規表現で「数字」に変換し、RMeCabで形態素解析をし、LDATS向けの形式のデータを作成していきます。
途中で、日本語の文字化け問題を回避するためにGoogle翻訳を使って単語名を置き換えています。
1時系列につき1文書となるようにデータを作っていく必要があるのですが、今回はボージョレ・ヌーボーのキャッチコピーなので最初から1時系列につき1文書となっているため都合が良いです。
データとソースコードはこちら

こちらは論文の図と同じものだとドキュメントの説明にあったので、論文の説明を見る限り、表すものとしては以下のようです。

  • 一番上の積み上げグラフはトピックごとの単語の割合を表しています。
  • 二番目の折れ線グラフはLDAによって推定されたトピックの時系列推移です。
  • 三番目のヒストグラムは二番目の時系列における変化点を集計したものです。
  • 四番目の折れ線グラフはモデルが推定したトピック割合の変化点の前後での推移です。

今回の図では文字が潰れていて見にくいですが、

  • トピック1はボキャブラリーが比較的リッチなコメント(「フルーティー」「フレグランス」「複雑」)
  • トピック2は数字を用いたコメント(「何年に一度の!」みたいな)
  • トピック3はボキャブラリーが貧相なコメント(「すごい!」みたいな)

のようです。
二番目の折れ線グラフを見る限り、周期的に数字を用いたコメントが現れているように思われます。四番目の折れ線グラフの変化点を見る限り、近年は数字を用いたコメントが相対的に減ってきて、リッチなボキャブラリーになってきているようです。

おわりに

時系列トピックモデルをカジュアルに試せる面白そうなパッケージだなと思い、LDATSパッケージを触ってみましたが、そもそもBoWなどを想定して作られているパッケージではないので、単語数が多いような分析ではそもそも可視化ができず使いにくいだろうなと思いました。マーケティングなどでユーザーのセグメントの推移を分析したい場合などにちょうど良いのだろうと思われます。

参考情報

[1] Long‐term community change through multiple rapid transitions in a desert rodent community
[2] Latent Dirichlet Allocation coupled with Bayesian Time Series analyses
[3] Package ‘LDATS’

RStudioをdockerで使える、Rockerを触ってみた

はじめに

dockerは会社のエンジニア向けの勉強会で紹介されていて、存在は知ってはいたが使っていませんでした。先日参加したTokyo.RでRockerというものが紹介されており、よし使ってみようと思うに至りました。
今回は遅ればせながら、Rockerを使ってみて、何かを回して、保存するという一連のチュートリアルをやってみます。

用語のざっくり理解

専門家の方に怒られたらあれですが、私の理解はこんな感じです。

  • Dockerfile
    • OS、言語、ライブラリ、バージョンなどをレシピの食材っぽい感じでテキスト形式で書かれたもの
  • イメージ
    • コンテナの元。Dockerfileをもとにdockerで生成できる。あるいはdocker hubとか言うのサイトから入手できる。ここからコンテナ(RStudioとかを回したい環境のこと)を起動できる。
  • コンテナ
    • RStudioとかを起動させておく実行環境。
  • 永続化
    • ファイルの内容が消えないように残すこと。立ち上げたコンテナを永続化しないと都度データやインストールしたパッケージが消える。

Rockerはdocker hubで提供されているイメージです。RockerのGitHubにはDockerfileも当然あげられていました。自分であれやこれや設定を追加するにはこれを参考にいじることになるのだと思います。

とにかくやってみる

まずは、dockerをインストールします。

Install Docker Desktop on Mac
docker自体はだいぶ前にインストールしたので、入れ方を忘れていますが、ここを見ればわかるはず。

ターミナルで以下を実行

ブラウザなどで以下を表示

DALEXのサンプルコードを回してみる。DALEXはTokyo.Rでも紹介されている機械学習の解釈可能性に関する手法をあらかた実践できるパッケージです。
moDel Agnostic Language for Exploration and eXplanation

以下の図はDALEXなどを使ってGLMで推定した結果の各特徴量の影響度をプロットしたものです。

ターミナルでコンテナのIDを確認する。

分析結果やソースコードやデータを保存したいので、コンテナの永続化を行います。

ここではユーザー名をperio、イメージ名をrocker_hadleyにしています。

永続化したコンテナのイメージがきちんと存在するか確認してみます。

ありますね。

さて、先程のコンテナを止めてみましょう。コンテナIDを指定するだけです。

永続化したコンテナを読み込んでRstudioを起動してみます。

実は、誤って最新の状態で保存せずに再度立ち上げてしまいましたwDALEXのサンプルコードを回すために色々パッケージインストールしたので、それが消えてしまった。残念。
元のimageにはないRDataがあるので、一応、永続化はできているので安心です。うむ、次からはきちんと保存をしよう。

これでどのPCでも同じ条件で分析していけるので、良さそうです。
ただ、一部で入らないパッケージ(sfパッケージとか)があったので、より良いDockerfileを探し求めたいとも思いますね。

u_riboさんがsfパッケージが動くdocker imageを教えてくれました。これまで、一つの環境にギュッとパッケージを詰め込んでやってきたんですが、目的に応じてdocker imageを選ぶのが良いのでしょう。

せっかくなので先日作成した事故物件分析用のコードを回してみましょう。

mapviewもインストールできて、問題なく使えました。u_riboさんありがとう。

参考情報

[1] The Rocker Project Docker Containers for the R Environment
[2] The Rocker Images: choosing a container
[3] rocker/tidyverse
[4] How to save data
[5] DALEXverse and fraud detection

2019年に読んだデータ分析系の本の振り返り(21+1冊)

はじめに

2020年、あけましておめでとうございます。年末に自分自身を振り返ろうと思ったのですが、結局データ分析と勉強しかしていないわけで、書籍を振り返ろうと思うに至りました。私の知り合いのデータサイエンティストはだいたい全冊持っているであろうと思われますが、良い本だと思うので思い出していただければ幸いです。

1.『ベイズモデリングの世界』(岩波書店)

基本的に階層ベイズモデルを使って、個体ごとの異質性を考慮した分析手法が提案されています。前半はオムニバス形式で様々な先生がモデルの適用について執筆されており、後半では伊庭先生による階層ベイズモデルの講義になっています。途中でスタイン統計量による縮小推定の話があげられ、柔軟なモデリングのためには「階層化した方が少なくとも望ましい推定量が得られる」という数学的証明を捨てることもやむを得ないと書かれています。

2.『トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ) 』(コロナ社)

この本はトピックモデルの教科書というよりも、ベイズ推定の教科書という側面が強い印象があります。途中で出てくる数式は流し読みするのは難しく、最低2冊以上のノートが別途必要になると思います。一度でもLDAのパラメータを導出してみたいという方には良い教科書だと思います。疑似コードが提供されているので、それをもとにRやPythonでコーディングしていけば、一番シンプルなLDAが非常に短い行で実行できてしまうことに驚かれるかもしれません。人間が手を動かして推定アルゴリズムを導出しているからこそ、短いコードで済むということを実感できるはずです。

3.『構造的因果モデルの基礎』(共立出版)

グラフィカルなアプローチで因果推論を扱っている書籍です。Judea Pearl流の因果推論アプローチについて記すことを目的に書かれています。基礎と書かれていますが決して簡単ではありません。ただ、扱われる数学のレベルとしては確率と線形代数がわかれば大丈夫だと思われます。余談ではありますが、1章の相関関係と因果関係の事例紹介で「おむつとビールの話」が都市伝説ではなくきちんと記事としてWall Street Journalという雑誌に掲載されていたことが明らかにされています。

4.『現場で使える!PyTorch開発入門 深層学習モデルの作成とアプリケーションへの実装 (AI & TECHNOLOGY)』(翔泳社)

PyTorchを触ったことがないが、深層学習の手法について知っている層を対象とした本です。6章まではGoogleのColabで動かせるのでGoogleに課金することなく深層学習による回帰、CNN、GAN、RNN、Encoder-Decoderモデル、ニューラル行列因子分解をPyTorchで試すことができます。写経したものはこちら。転移学習や高解像度化や画像生成、文章のクラス分類、文書生成、機械翻訳などもできるので、PyTorchでこれくらいの量をコーディングしたらこれくらいのことができるのかという学びや、他の人の書いたPyTorchコードを読みやすくなるなどの便益は十分にあると思いました。

5.『作ってわかる! アンサンブル学習アルゴリズム入門』(シーアンドアール研究所)

会社で行っているPythonもくもく会用に買った本で、scikit-learnを使わずに機械学習のアルゴリズム(アンサンブル系)をコーディングするための本です。pythonのコードについて逐次、細かい解説が行われているわけではないので、1行1行自分でコメントを加えながら写経をしていけば力が付くという本かなと思われます。sklearnはそれはそれで素晴らしいですが、こういう本でフルスクラッチで修行できるのはいいですね。

6.『数理統計学―基礎から学ぶデータ解析』(内田老鶴圃)

統計検定1級を合格された方のブログで紹介されていた教科書です。理系の大学生レベルの数学知識があれば、数理統計学の基礎を学べると思います。中心極限定理の証明や、様々な分布の期待値や分散、様々な分布の性質について数式を用いてしっかり理解することができます。数式もほどよく端折られているので、無論ですがノートが数冊必要になります。各章毎にある練習問題も解くことで力が付くと思います。日本の大学の授業の教科書がこれだったらジェノサイド(再履修者の大量発生)が起きるんだろうなと思ってしまった。

7.『44の例題で学ぶ統計的検定と推定の解き方』(オーム社)

統計の検定に関してだけ扱った珍しい本です。第3部までは統計学の普通の教科書ですが、それ以降であらゆる検定の例題が44件も載せられています。パラメトリックな検定から、ノンパラメトリックな検定まで幅広く扱われています。一番気にいっているのは仮説検定法の分類の表です。これさえあれば、どのデータに対してどの検定を行えばいいかが一目瞭然です。

8.『わけがわかる機械学習 ── 現実の問題を解くために、しくみを理解する』(技術評論社)

機械学習の原理を手早く数式を交えて学べる本です。かゆいところに手が届いていると言うか、既出の教科書では捨象されがちな、条件付き確率における2変数以上の条件づけでの表現に紙面を割いていたりしてくれるのが嬉しいです。ある程度数学の話はわかるが、だいぶ忘れているビジネスパーソンには大変にありがたいコンテンツと言えると思います。ベイズ線形回帰に関しても行列を用いた、わかりやすい導出方法が紹介されています。またコラムで紹介されている、測度論にどう向き合えばいいかの著者の見解は参考になります。

9.『Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)

R言語とstanを用いてベイズ統計学を入門レベルから学べる本です。各トピックごとにそれなりの紙面が割かれています。例題も豊富にあるので、線形回帰・MCMC・情報量基準・階層ベイズモデルまで、ベイズ統計学を基礎から応用までしっかりと学べると思います。youtubeで著者の講義も配信されているので、留学気分を味わえます。

10.『scikit-learnとTensorFlowによる実践機械学習』(オライリージャパン)

2019年に日本で開かれたML SummitでTFの開発者がおすすめしていた教科書です。前半部分で機械学習の入門から応用までをわかりやすい説明で学ぶことができます。数式は少ないですが、図とソースコード(Python)がちりばめられており、手を動かして理解を進めることができます。後半部分はTensorFlowを用いた深層学習の基礎を同様に手を動かして学ぶことができます。ただ、TFのバージョンも変わってきているので前半の説明をアテにして読むのも良いと思います。

11.『AIアルゴリズムマーケティング 自動化のための機械学習/経済モデル、ベストプラクティス、アーキテクチャ (impress top gear)

マーケティングへのデータサイエンスの適用に関する珍しい書籍です。ソースコードはついていないですが、業務で使う際のアイデアが手に入ることもあります。一般的な回帰、生存時間分析、オークション、アトリビューション分析、アップリフトモデリング以外にも、情報検索やレコメンデーションやトピックモデルなどマーケティングながら学際的なトピックも扱われています。レコメンドなどで使われる、ランク学習に関して詳しく書かれた書籍をあまり知らないので、この本はその点においてもありがたい本でもあります。

12.『入門 統計的因果推論』(朝倉書店)

ほぼ全ての章でグラフィカルなアプローチで因果推論を扱っています。例題も豊富なので、一つ一つ丁寧にやれば理解が捗ります。おそらく、例題の多さを含め一番丁寧にd分離性、do演算子、バックドア基準、フロントドア基準に関する説明をしてくれている本なのかなと思いました。グラフでの因果推論に関して初めての人でも、確率さえ知っていれば読み進めることができるはずです。また、途中で操作変数法の紹介もされ、経済学出身者としては読みやすい。ただ、傾向スコアのくだりや、DIDなどのくだりはあまり出てきません。あと、やってないですが章末の練習問題に対するSolution Manualが提供されているようです。

13.『実践 ベイズモデリング -解析技法と認知モデル-』(朝倉書店)

ベイズモデリングを様々な事例に適用する方法がオムニバス形式で記された本です。ワイブル分布や異質性を考慮した二項分布、無制限複数選択形式のアンケートデータに対する手法、トピックモデル、項目反応理論などが扱われています。マーケティングの実務で使える事例が多いように感じました。こちらはサポートサイトでRコードとstanコードが提供されています。あと、appendixにあるプレート表現の見方も参考になります。

14.『機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

機械学習などで用いるベイズ推論を扱った教科書です。入門とありますが、入門者は書かれた数式をそのまま見ていても頭に入らないのではないでしょうか。手を動かしてなんぼの本だと思います。ノート2冊は絶対に必要です。たぶん、数式の展開を丁寧に記すと倍以上の厚みの本になると思います。各々のモデルに関してグラフィカルモデルが記されているのや、サンプルコードとしてGitHubにJuliaで書かれたソースコードが提供されているのも良いです。

15.『その問題、数理モデルが解決します』(ベレ出版)

物語形式で、様々な問題に対して数理モデリングのアプローチが紹介されています。途中でマッチング理論やゲーム理論やオークションなども登場することから、経済学出身者も喜ぶ内容かもしれません。社会人になってからナッシュ均衡という言葉が書かれた本は中々出会って来なかった。

16.『ヤバい予測学 ― 「何を買うか」から「いつ死ぬか」まであなたの行動はすべて読まれている』(CCCメディアハウス)

2013年と結構古い本ですが、データ分析を様々な事象に対して適用した事例紹介本です。アップリフトモデリングへの言及もあり、こういったものに関して日本は何年も遅れてブームが来るんだなという実感を与えてくれた本でもありました。appendixに分析事例が147個ほどあげられているのも参考になります。

17.『たのしいベイズモデリング2: 事例で拓く研究のフロンティア』(北大路書房)

主にstanを用いたベイズモデリングによる分析事例が1と2で38本もオムニバス形式で載っています。ほとんどの事例で階層ベイズモデルが扱われています。2では若干マーケティングに近い内容の題材も扱われ、データサイエンティストの人にも嬉しい内容かもしれません。もちろんデータとstanとRのコードがサポートサイトで提供されています。

18.『カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― (統計学One Point 2)』(共立出版)

状態空間モデルで時系列予測を行うための手法が記されている本です。RのKFASパッケージが全面に渡って扱われています。トレンドを考慮したり、カレンダー効果を追加したり、共変量を追加したりなど様々なアプローチが紹介されコードも伴っているわけですから、業務でも抜群に役に立ちました。

19.『機械学習のエッセンス -実装しながら学ぶPython,数学,アルゴリズム- (Machine Learning)』(SBクリエイティブ)

自分のいる会社で最低限の数学がわかると思われる若いメンバーに買ってもらうように言っている本です。微積分・線形代数だけでなく、カルシュ・キューン・タッカー条件(最適化数学)に関しても扱ってくれているので、ここで出てくる数学がわかれば大体の論文に立ち向かえると思います。さらに、Pythonの基礎もこれで学ぶことができるので一石二鳥な素敵な本ですね。また、最後の方でスクラッチでアルゴリズムを書くパートがあり、こちらも勉強になります。

20.『機械学習のための特徴量エンジニアリング ―その原理とPythonによる実践 (オライリー・ジャパン)』(オライリー・ジャパン)

機械学習における前処理の指針を与えてくれる本です。Pythonのコードが提供されています。例えばですが、「テストデータにだけある、新しい単語は取り除いてしまえばいい」などの細かいアドバイスが何気に嬉しいです。「Effectコーディング」「特徴量ハッシング」「ビンカウンティング」「バックオフ」「leakage-proof統計量」などは読むまで知らないところだったので勉強になりました。

21.『データサイエンスのための統計学入門 ―予測、分類、統計モデリング、統計的機械学習とRプログラミング』(オライリージャパン)

データ分析の仕事をする上で最低限必要な知識を幅広く抑えることができる本です。数式は少なく、ところどころ出てくるコードはR言語です。参考文献などがブログだったりするため厳密さがめちゃあるわけではないですが、業務で使う分には問題ないと思います。分類問題において、AUCなどの評価指標だけでなく、予測値自体の探索的分析のすすめなどが書かれており、参考になりました。また、特徴量エンジンとしてのk-NN法の話も面白いと思いました。

[+α]『プログラマのためのGoogle Cloud Platform入門 サービスの全体像からクラウドネイティブアプリケーション構築まで』(翔泳社)

Google Cloud Platformを初めて触るデータ分析者にはちょうど良い本です。説明もわかりやすいので、いきなりアカウントを作ってドキュメントを解読するよりかは戸惑いは減るはずです。この本を土台に、GCS・GCEを駆使してML系のAPIを呼び出して使うなどの最低限の操作は私でもできるようになりました。GCPの画面や機能もどんどん変わっていくので書籍を買ってもアレなんですが、歴史的な背景も若干記述されているので、それはそれで勉強になります。ただ、エンジニアにこの本を買うべきか聞いた際にネガティブな意見があったのですが、たぶん現役プログラマからすると簡単過ぎるからなんだろうなと思います。

終わりに

2019年もぼちぼち勉強できましたが、2020年もこれまで同様にノートとペンを大事にする勉強を続けていき、コーディングも分析ももっともっと数をこなして会社や社会に求められるようなデータ分析官を目指していこうと思います。あぁ、英会話などの勉強をする時間を作るのが難しい。

R advent calendar 2019 RSelenium、jpmesh、sfパッケージで東京23区の事故物件を分析してみよう!

はじめに

今回で3回目となるR advent Calendarですが、前回は「一発屋芸人の検索トレンドのデータ」を扱い、前々回は「ポケモンのデータ」を扱いました。今回は人の命に関わるようなデータを扱ってみたいと思い、某サイトから東京都の23区内における事故物件の住所と詳細を集めてきました。どのようなエリアが事故が起きやすいのかの分析を行います。(以下では、事故物件をAP(Accident Property)と呼びます。)

分析工程

・データの収集
・データの整形
・可視化
・分析

データの収集

APに関する情報を某サイトより集める必要があります。そこで必要なライブラリとしては、RSeleniumやtwitteRがあげられます。
twitteRが必要な理由は、APに関するサイトにAPの一覧ページがなく、公式アカウントがAPに関するページのリンクを投稿しているところにあります。ただ、私が以前使っていた頃とはTwitterAPIの仕様が変わり、3ヶ月よりも前の情報にアクセスできなくなっていました。そのため、今後のデータに関してはTwitterAPIでいいのですが、過去のものに関しては別アプローチが必要となります。
また、APに関するサイトはJavaScriptで地図が表示されているだけなので、RSeleniumを使って地図をクリックさせ、表示された情報をスクレイピングするなどの処理が必要となります。
当初の想定ではTwitterのデータ収集、リンク先の情報をRSeleniumでスクレイピングするだけの簡単な仕事だと思っていたのですが、過去のデータにアクセスできないので、地図上で一つ一つ見つけていくためにRSeleniumだけで頑張ることにしました。(私の過去のアドカレ史上、一番面倒なデータとなりました。)

誰もやらないと思うのですが、一応手順を記しておきます。

RSeleniumだけでMapからAPの情報を抽出するための手順
1.都内の住所一覧を収集
2.検索窓に住所を入力
3.検索結果一覧の上位5件をクリック
4.一度地図を引くことでAPを広い範囲で捉えれるようにする
5.APのマークの要素を取得し、1件ずつクリックし、表示されたAPの情報をデータフレームに格納する

こちらが取得できたデータです。

RSeleniumでAPの
・住所
・発生時期(フリーテキスト)
・AP詳細
・投稿日
を集めることができるので、その住所データに対して、Yahoo!のジオコードに関するAPIを利用します。(利用申請が必要なはずです。4年前くらいに申請していたのでそのまま使えました。)
Yahoo!のAPIを使えば、住所から緯度経度の情報を取得することができます。

APの緯度経度がわかれば、jpmeshパッケージを用いて1kmメッシュやら10kmメッシュやらのメッシュデータに変換することができます。
jpmeshを用いてメッシュデータに変換し、メッシュ単位でAPの発生件数を集計します。

データ収集用のソースコードは思いのほか長くなってしまったので、GitHubにあげておきました。
https://github.com/KamonohashiPerry/r_advent_calendar_2019

ここで再度、手順を整理しておきましょう。

Twitterに出てきたものだけを取得(直近3ヶ月〜)する場合、
run_tweet_collect.RでTweetを収集

run_selenium.RでAPの情報をスクレイピング

run_map_api.Rで住所から緯度経度の取得

making_mesh_data_and_download_other_data.Rで緯度経度からメッシュデータへの変換、その他の人口データや地価データと接続をします。

直近3ヶ月以前のものを取得する場合、
run_selenium_from_map.Rで地図上から直接APの情報を取得する

making_mesh_data_and_download_other_data.Rで緯度経度からメッシュデータへの変換、その他の人口データや地価データと接続をします。

データの整形

APの1kmメッシュデータを手に入れたら、kokudosuuchiパッケージを使って国土地理院の収集したデータをつなぎこみます。手順としては、以下のとおりです。

まずは推計人口というそのエリアの人口の予測値です。今回は2010年のものを抽出しました。こちらは1kmメッシュのデータなので、変換することなく使えて都合が良いです。

続いて、2015年の公示地価を抽出しました。

こちらはメッシュデータではないので、緯度経度の情報から1kmメッシュのデータに変換する必要があります。後で行います。

単純集計・可視化

今回のデータセットのデータ数は3919件です。本当は7000件以上はあると思われますが、マップから取ってくるという勝手上、なかなか全てを取り切ることができませんでした。

まずは、1kmメッシュごとのAP発生件数のヒストグラムです。

1kmメッシュにおける人口のヒストグラムです。

公示地価のヒストグラムです。

1kmメッシュにおける人口あたりのAP件数のヒストグラムです。

分析

ここでは、色々な軸でAPのデータに向き合ってみようと思います。

APの発生件数の集計

人口が多いところがAPの発生件数が多いところだと思われますが、とりあえず確認します。

世田谷区は最も人口が多いことから、AP発生件数では一番となっています。続いて、歌舞伎町などがある新宿が来ています。しかしながら、人口に占めるAP発生件数で言うと、港区がかなり高く出ているのがわかります。

人口あたりのAP件数

ここでは、メッシュデータをsfパッケージ用のオブジェクトに変換して、1kmにおける人口あたりのAP発生割合を可視化しています。

こちらのmapviewパッケージで作ったマップはインタラクティブにいじることができます。ぜひ関心のあるエリアでいじってみてください。

1kmメッシュ人口あたりのAP発生件数(×100)の可視化

1kmメッシュでのAP発生件数の可視化

比率ベースで、色の明るいメッシュのところを見ると、港区、中央区、新宿区、渋谷区などがAPが発生しやすいようです。件数ベースで言うと新宿が一番多いですね。
一番色が明るい港区はてっきり六本木ではないかと思ったのですが、新橋から日比谷にかけたエリアでした。会社員による自○が多いようです。恐ろしいものです。

APの名前の集計

APの名前を集計してみます。これは別にこの名前だからAPになりやすいというわけではなく、単純に数が多いだけの可能性がありますし、実際にそうだろうと思われます。AP発生率を知るには、APではないものも含めた全物件名に占めるAP発生件名を手に入れないといけませんが、全物件名を収集するのが難しいことから単純に頻度の集計となります。今回は、wordcloud2パッケージを使って、ワードクラウドにしてみます。文字が大きいと頻度が高いものとなります。

ハイツ、荘、コーポ、マンション、号棟、アパート、ハウスなどが多く出現しているようです。ただ、物件の名前としても頻度が高いとも考えられますね。

地価と人口あたりのAP発生件数の関係

ここでは地価のデータとAP発生の関係性について見てみます。

地価の階級値(10個のパーセンタイルに分割)を横軸に、縦軸に人口あたりAP発生数をおくと、地価が上がるに従い人口あたりAP発生数が高まる傾向があります。これは、人口密度が高く地価の高いところではAPが発生しやすいということを示しているのではないでしょうか。人口密度が高いと地価があがる、人口密度が高いと治安が悪くなるという可能性が考えられます。

APの詳細の集計

ここではAPになってしまった詳細の内容について先ほどと同様に形態素解析を行いワードクラウドにしてみます。

どうやら孤独死が多いようです。高齢者の人口構成比が関係しているのだろうと思われます。

APの発生時期に関するテキストマイニング

ここでは、発生時期に含まれる四桁の数字を集計して、何年くらいのAPが多いのかをざっくりと掴みます。

どうやら昔のデータはあまり登録されていないようです。記憶が確かではないかもしれませんし、古すぎるものは消されているのかもしれませんね。あのサイトはユーザー生成コンテンツ(UGC)なので、投稿する人はそこまで昔のことをわざわざ投稿しないのかもしれないですね。

APの詳細に関するテキストマイニング

ここではトピック数10として、topicmodelsパッケージを使いLDAを行います。

なかなか恐ろしいキーワードが多いですが、なんとなくですがうまく分類されているのではないかと思われます。

トピック1は男性の不幸
トピック2は不動産屋に言われた告知事項
トピック3は孤独死
トピック4は病死
トピック5は火災・転落・事故
トピック6は事故のあった建物に関する記載
トピック7は腐乱した事例
トピック8は建物に関して不明であることの記載
トピック9は心理的瑕疵あり
トピック10は自○

となっているように思われます。まさかこのようなデータにトピックモデルを使うことになるとは。

おわりに

今回はR言語のみを用いて、APに関するデータを収集し、地図にプロットしたり他のメッシュデータとつなぎ合わせて分析をするなどしました。APが発生しやすいエリア、APと地価との関係、APのテキストマイニングなど興味深い結果が得られたと思います。
一つ残念なのは、時系列情報がフリーテキストなので、APがどのエリアでどの頻度で発生していくのかの分析のコストが高く、今回は時系列情報を用いた分析にチャレンジできませんでした。
今後はタクシーの需要推定の分析で行われているように、メッシュ単位でのAP発生確率の推定などを機械学習で行えると面白いなと思います。どなたか一緒にアノテーションしましょう!

それでは、どうか良い年末をお過ごし下さい!
メリークリスマス!

参考情報

数多くの方々の記事を見てどうにか仕上げることができました。感謝します。

[1]【追記あり】sfパッケージでシェープファイルを読み込んでmapviewパッケージで可視化するまで
[2]How to use mesh cord in R
[3]Rを使ってワードクラウドを作ってみました
[4]国土数値情報ダウンロードサービスWeb APIからデータを取得するためのRパッケージです
[5]東京の地価公示データを眺める
[6]Chapter 1 Introduction to spatial data in R
[7][翻訳] RSelenium vignette: RSeleniumの基本
[8]RからYahoo!のジオコーディングを利用する方法
[9]EMBEDDING A LEAFLET MAP ON WORDPRESS
[10]mapview advanced controls
[11]RSeleniumでChromeからファイルをダウンロードするディレクトリを指定する方法
[12]Selenium Serverが立ち上がらないときはportが被っているかも!?
[13]brew install selenium-server-standalone
[14]ナウでヤングなRの環境変数管理方法
[15]タクシードライバー向け需要予測について
[16]LDA with topicmodels package for R, how do I get the topic probability for each term?
[17]dplyr — 高速data.frame処理

[数理統計学]連続型確率分布の期待値と分散の導出まとめ

はじめに

前回は、離散型確率分布に関するよくある分布とその平均と分散の導出についてひたすら記しました。今回は連続型の確率分布に関して同様に記していこうと思います。

※PC版でないと数式が切れてしまうので、SP版での閲覧はおすすめしません。私、SKUEのスマホだけは表示崩れがないようにチェックしております。通勤電車で定期的に読むので。

今回登場する連続分布

『数理統計学―基礎から学ぶデータ解析』という本に登場するものを式を補いながら紹介します。

1.一様分布

離散型の一様分布に対応しており、

$$ f(x; \alpha, \beta) = \frac{1}{\beta – \alpha}, \alpha \leq x \leq \beta $$
で定義されます。

ここで、ある区間(a,b)が区間(α,β)に含まれるとした場合、Xが一様分布に従う際の確率は、

$$ P(X \in (a,b)) = \int_a^b \frac{1}{\beta – \alpha} dx = \frac{b-a}{\beta – \alpha} $$
となり、確率は区間の長さ(b-a)に比例し、その区間が(α,β)のどこにあってもよい。つまり、一様にどこでも一定の確率となることを示している。

また、一様分布から生成される変数は、ある確率分布に従うような乱数を、確率分布関数の逆数から生成させることができる。そのような乱数生成方法を逆関数法と呼ぶ。

逆関数法
[0,1]上の一様分布に従う確率変数をUとする。このとき、目的の確率分布の関数の逆関数に一様乱数Uを代入したものは、目的の確率分布に従う。

Uが一様分布に従うことを式で示すと、
$$ P(U \leq u) = u, 0 \leq u \leq 1$$
ここで、
$$ U \leq u \Leftrightarrow F^{-1}(U) \leq F^{-1}(u) $$
であるから、
$$ P( F^{-1}(U) \leq F^{-1}(u) ) = u $$
となる。
さらに、\( F^{-1}(u) = x \)とおくと、
$$ P( F^{-1}(U) \leq x ) = F(x) $$
これは確率変数、\( F^{-1}(U) \)が累積分布関数が、\( F(x) \)であるような確率分布に従うことを表している。

一様分布の平均

$$ E(X) = \int_a^b xf(x) dx $$
$$ = \int_a^b x \frac{1}{b-a} dx$$
$$ = \left[ \frac{x^2}{2} \frac{1}{b-a} \right]^b_a $$
$$ = \frac{b^2 – a^2}{2(b-a)} $$
$$ = \frac{(b-a)(b+a)}{2(b-a)} $$
$$ = \frac{b+a}{2} $$

一様分布の分散

$$ E(X^2) = \int_a^b x^2f(x) dx $$
$$ = \int_a^b x^2 \frac{1}{b-a} dx $$
$$ = \left[ \frac{x^3}{3} \frac{1}{b-a} \right]^b_a $$
$$ = \frac{(b-a)(b^2 + ab + a^2)}{3(b-a)} $$
$$ = \frac{(b^2 + ab + a^2)}{3} $$

$$V(X)=E(X^2)-[ E(X) ]^2$$
$$ = \frac{(b^2 + ab + a^2)}{3} – \left( \frac{b+a}{2} \right)^2$$
$$ = \frac{(b^2 – 2ab + a^2)}{12} $$
$$ = \frac{(b – a)^2}{12} $$

2.正規分布

正規分布は以下のような式で定義される確率分布のことを指す。

$$ f(x; \mu, \sigma^2) = \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(x-\mu)^2}{2\sigma^2} } $$
$$ -\infty < x < \infty , -\infty < \mu < \infty, \sigma > 0 $$

上の式で、μ=0、σ=1の場合は標準正規分布と呼ぶ。

実際に、正規分布に従う変数Xを標準化したもの\( Z = \frac{X-\mu}{\sigma} \)が標準正規分布に従うことを以下で示す。
$$ P(Z \leq z) = P(X \leq \mu + \sigma z) $$
$$ = \int_{-\infty}^{\mu + \sigma z} \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(x-\mu)^2}{2\sigma^2} } dx $$
ここで\( y = \frac{x-\mu}{\sigma} \)とすると、
$$ = \int_{-\infty}^{z} \frac{1}{ \sqrt{2 \pi} } e^{ – \frac{y^2}{2} }dy $$
となり、標準正規分布に従っていることがわかる。最後の積分の上限はもともとのxに対して、μを引き、σで割るという計算をすることで得られる。

また、正規分布の平均値や分散を計算する上で避けられない、ガウス積分について記しておく。

まず、
$$ I = \int_{-\infty}^{\infty} e^{-\frac{x^2}{2}} dx $$
とおき、xとyについての積分で

$$ I^2 = \left[ \int_{-\infty}^{\infty} e^{-\frac{x^2}{2}} dx \right] \left[ \int_{-\infty}^{\infty} e^{-\frac{y^2}{2}} dy \right] $$

$$ = \int_{-\infty}^{\infty} \!\!\! \int_{-\infty}^{\infty} e^{-\frac{x^2+y^2}{2}}dxdy $$
となるようにする。
ここで、重積分の変数変換を行う。\( x = r \cos \theta \)、\( y = r \sin \theta \)とし、dxdyをdθdrに変換するためにヤコビアンを用いて、
$$ = \int_{-\infty}^{\infty} \!\!\! \int_{-\infty}^{\infty} e^{-\frac{r^2}{2}} \left| \begin{array}{ccc} \frac{\partial x}{\partial r} & \frac{\partial y}{\partial r} \\ \frac{\partial x}{\partial \theta} & \frac{\partial y}{\partial \theta} \end{array} \right| d\theta dr $$
と置き換える。

$$ = \int_{-\infty}^{\infty} \!\!\! \int_{-\infty}^{\infty} e^{-\frac{r^2}{2}} \left| \begin{array}{ccc} \cos \theta & \sin \theta \\ -r \sin \theta & r \cos \theta \end{array} \right| d\theta dr $$

$$ = \int_{-\infty}^{2\pi} \!\!\! \int_{0}^{\infty} e^{-\frac{r^2}{2}} r d\theta dr $$

$$ = \left[ \theta \int_{0}^{\infty} e^{-\frac{r^2}{2}} r dr \right]^{2\pi}_0$$
$$ = 2\pi \int_{0}^{\infty} e^{-\frac{r^2}{2}} r dr $$
$$ = 2\pi \left[ – e^{ -\frac{r^2}{2} } \right]^{\infty}_0$$
$$ = 2\pi $$
よって、
$$ I = \sqrt{2\pi} $$
となる。

さっそく、ガウス積分を使って、正規分布の積分が1になることを示す。

$$ \int_{-\infty}^{\infty} f(x) dx $$
$$ = \int_{-\infty}^{\infty} \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(x-\mu)^2}{2\sigma^2} } dx $$
ここで、\( z = \frac{x-\mu}{\sigma}\)とおき、\( \frac{dx}{dz}=\sigma\)から、
$$ = \int_{-\infty}^{\infty} \frac{1}{ \sqrt{2 \pi} } e^{ – \frac{z^2}{2} } dz $$
となる。
ガウス積分の公式(\( \int_{-\infty}^{\infty} e^{-\frac{z^2}{2}} dz = \sqrt{2\pi} \) )より、

$$ = \frac{1}{\sqrt{2 \pi}} \times \sqrt{2\pi} = 1 $$

正規分布の平均

$$ E(X) = \int_{-\infty}^{\infty} x \times \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(x-\mu)^2}{2\sigma^2} } dx $$
ここで、\( z = \frac{x-\mu}{\sigma}\)とおき、\( \frac{dx}{dz}=\sigma\)から、
$$ = \int_{-\infty}^{\infty} (\mu + \sigma z) \times \frac{1}{ \sqrt{2 \pi} } e^{ – \frac{z^2}{2} } dz $$
第一項目は定積分が1になることから、第二項目は奇関数であることから定積分が0になるため、
$$ = \mu $$
となる。

正規分布の分散

$$ V(X) = E(X^2) – \left[ E(X) \right]^2 $$
$$ = \int_{-\infty}^{\infty} x^2 \times \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(x-\mu)^2}{2\sigma^2} } dx – \mu^2 $$
ここで、\( z = \frac{x-\mu}{\sigma}\)とおき、\( x = \sigma z + \mu \)、\( \frac{dx}{dz}=\sigma\)から、
$$ = \int_{-\infty}^{\infty} (\sigma z + \mu)^2 \times \frac{1}{ \sqrt{2 \pi}} e^{ – \frac{z^2}{2} } dz – \mu^2 $$

$$ = \int_{-\infty}^{\infty} (\sigma^2 z^2 + 2\sigma \mu z + \mu^2) \times \frac{1}{ \sqrt{2 \pi}} e^{ – \frac{z^2}{2} } dz – \mu^2 $$

ここで、第二項目は奇関数であることから定積分が0になるため、第三項目は定積分が1になることから、第三項目と第四項目が消し合うため、
$$ = \int_{-\infty}^{\infty} \sigma^2 z^2 \frac{1}{ \sqrt{2 \pi}} e^{ – \frac{z^2}{2} } dz $$
となる。
部分積分の公式より、

$$ = \frac{\sigma^2}{\sqrt{2 \pi}} \left( \left[ \frac{z^2}{-z} e^{ – \frac{z^2}{2} } \right]^\infty_-\infty – \int_{-\infty}^{\infty} 2z \times \frac{1}{-z} e^{ – \frac{z^2}{2}} \right) $$

第一項目はガウス積分に他ならないため、
$$ = \frac{\sigma^2}{\sqrt{2 \pi}} ( – \sqrt{2 \pi} + 2 \sqrt{2 \pi} ) $$
となり、
$$ = \sigma^2 $$
となる。

3.対数正規分布

正の値をとる確率変数Xを対数変換したものが正規分布に従うとした分布を対数正規分布と呼ぶ。

$$ f(x; \mu, \sigma^2) = \frac{1}{ \sqrt{2 \pi} \sigma x } e^{ – \frac{( \log x-\mu)^2}{2\sigma^2} } , x > 0$$

対数正規分布は以下のように正規分布から導出することができる。

\( f_x(x) = \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(x-\mu)^2}{2\sigma^2} } \)とし、\( x = g(y) = e^y \)とおき、以下の分布関数について以下のように変形する。

$$ F_x (x) = \int_{-\infty}^{x} f_x(t)dt $$
$$ = \int_{-\infty}^{x} f_y(g^{-1}(x)) \frac{dt}{dx} dx $$
$$ = \int_{-\infty}^{x} f_y(g^{-1}(x)) \frac{d g^{-1}(x) }{dx} dx $$

ここで両辺をxで微分すると、

$$ f_x(x) = f_y(g^{-1}x) \frac{d g^{-1}(x) }{dx} $$
$$ = f_y( \log x) \frac{1}{x} $$
$$ = \frac{1}{ \sqrt{2 \pi} \sigma x } e^{ – \frac{( \log x-\mu)^2}{2\sigma^2} } $$

対数正規分布の平均

$$ E(X) = E(e^y) $$
$$ = \int_{-\infty}^{\infty} e^y g(y) dy $$
$$ = \int_{-\infty}^{\infty} e^y \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(y-\mu)^2}{2\sigma^2} } dy $$
$$ = \int_{-\infty}^{\infty} \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{( y – \mu – \sigma^2 )^2}{2\sigma^2} + \mu + \frac{\sigma^2}{2} } dy $$
ここで、\( y – \mu – \sigma^2 = t \)とすると、
$$ = e^{\mu + \frac{\sigma^2}{2}} \frac{1}{ \sqrt{2 \pi} \sigma } \int_{-\infty}^{\infty} e^{-\frac{t^2}{2\sigma^2} } dt $$
となり、定積分の値はガウス積分の公式より\( \sqrt{2\pi} \sigma \)なので、
$$ = e^{\mu + \frac{\sigma^2}{2}} $$

対数正規分布の分散

$$ E(X^2) = E(e^{2y}) $$
$$ = \int_{-\infty}^{\infty} e^{2y} g(y) dy $$
$$ = \int_{-\infty}^{\infty} e^{2y} \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(y-\mu)^2}{2\sigma^2} } dy $$
$$ = \int_{-\infty}^{\infty} \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{1}{2\sigma^2} \left( y – (\mu + 2\sigma^2 ) \right)^2 + 2\mu + 2\sigma^2 } dy $$
ここで、\( y – 2\sigma^2 = s \)とすると、
$$ = e^{2\mu + 2\sigma^2} \frac{1}{ \sqrt{2 \pi} \sigma } \int_{-\infty}^{\infty} e^{ – \frac{(s-\mu)^2}{2\sigma^2} } ds $$
ここで、\( s – \mu = t \)とおき、ガウス積分の公式を用いると、
$$ = e^{2\mu + 2\sigma^2} $$

$$ V(X) = E(X^2) – \left[ E(X) \right]^2 $$
$$ = e^{2\mu + 2\sigma^2} – e^{2\mu + \sigma^2} $$
$$ = e^{2\mu + \sigma^2} (e^{\sigma^2} – 1) $$

4.ガンマ分布

$$ f(x; \alpha, \beta) = \frac{\beta^{\alpha} }{ \Gamma (\alpha) } x^{\alpha – 1}e^{-\beta x}, \alpha > 0, \beta > 0, x >0 $$

と表される確率密度関数をもつ分布をガンマ分布と呼ぶ。
\( \alpha, \beta \)の2つのパラメータを持ち、様々な形になり、正の値をとる変数の確率分布のモデルとして使われる。

ガンマ分布はポアソン過程から導出することができる。その際は、ガンマ関数が階乗の一般化であるという性質を使う。

まずガンマ関数が階乗の一般化であることの証明を行う。

まず、ガンマ関数は以下のように定義される。
$$ \Gamma (\alpha) = \int_{0}^{\infty} e^{-t} t^{\alpha – 1} dt $$
\( \alpha = 1 \)の場合、
$$ \Gamma (1) = \int_{0}^{\infty} e^{-t} dt $$
$$ = \left[ -e^{-t} \right]^{\infty}_{0} = 1 $$

任意の正の整数nに対して、
$$ \Gamma (n ) = \int_{0}^{\infty} e^{-t} t^{ n – 1} dt $$
部分積分の公式より、

$$ = \left[ -e^{-t} t^{n-1} \right]^{\infty}_{0} + (n-1)\int_{0}^{\infty} e^{-t}t^{n-2}dt$$

第1項目は0に、第2項目はガンマ関数の定義から以下のようになる。
$$ = (n-1) \Gamma (n-1) $$
以上より、逐次的に代入することで、
$$ \Gamma (n+1) = n! \Gamma (1) = n! $$
となる。

以下、ポアソン過程からガンマ分布を導出する。
まず、分布関数
$$ F_n(t) = P( Y_n \le t ) = P(X_t \ge n) $$
を考える。tや\( Y_n \)はn回の事象が生じるまでの時間を表し、\( X_t \)は事象の回数を表している。
\( F_n \)はポアソン分布に従うため、
$$ = \sum _{x=n}^{\infty} \frac{e^{-\lambda t} (\lambda t)^x }{x!} $$
と表され、全事象からn-1回分までの確率を差し引くことによっても表される。
$$ = 1 – \sum _{x=0}^{n-1} \frac{e^{-\lambda t} (\lambda t)^x }{x!} $$

\(Y_n \)の密度関数は分布関数をtで微分することによって得られる。

$$ \frac{d F_n (t)}{dt} = \sum _{x=0}^{n-1} \frac{\lambda e^{-\lambda t} (\lambda t)^x }{x!} – \sum _{x=1}^{n-1} \frac{\lambda e^{-\lambda t} \lambda^x t^{x-1} }{(x-1)!} $$
最後の項以外はキャンセルアウトされるため、
$$ = \frac{\lambda^n t^{n-1}e^{-\lambda t}}{(n-1)!} $$
となり、先程示したガンマ関数が階乗の一般系であることの証明から、
$$ = \frac{\lambda^n}{\Gamma (n) } t^{n-1}e^{-\lambda t} $$
となり、これはガンマ分布である。つまり、ポアソン過程に従う確率分布を時間で微分することでガンマ分布は得られる。

また、ガンマ分布には再生性があり、モーメント母関数を用いたり、確率分布の畳込みを用いての証明などがある。

ガンマ分布の平均

$$ E(X) = \int_{0}^{\infty} x \frac{\beta^{\alpha} }{ \Gamma (\alpha) } x^{\alpha – 1}e^{-\beta x} dx $$
$$ = \int_{0}^{\infty} \frac{\beta^{\alpha} }{ \Gamma (\alpha) } x^{\alpha}e^{-\beta x} dx $$
部分積分の公式より、

$$ \frac{ \beta^{\alpha} }{ \Gamma ( \alpha )} \left[ \frac{ x^{\alpha} e^{-\beta x} }{ -\beta } \right]^{\infty}_0 – \frac{ \beta^{\alpha} }{ \Gamma ( \alpha )} \int_{0}^{\infty} \frac{ \alpha x^{\alpha – 1}e^{-\beta x}}{-\beta} dx $$

第1項目は0になるため、
$$ = \frac{\alpha}{\beta} \int_{0}^{\infty} \frac{\beta^{\alpha} x^{\alpha-1}e^{-\beta x}}{\Gamma (\alpha) } $$
となり、ガンマ分布の定積分は1であることから、
$$ = \frac{\alpha}{\beta} $$
となる。

ガンマ分布の分散

$$ E(X^2) = \int_{0}^{\infty} x^{2} \frac{\beta^{\alpha} }{ \Gamma (\alpha) } x^{\alpha – 1}e^{-\beta x} dx $$
$$ = \frac{\beta^{\alpha}}{\Gamma (\alpha) } \int_{0}^{\infty} x^{\alpha + 1}e^{-\beta x} dx $$
部分積分の公式より、

$$ = \frac{\beta^{\alpha}}{\Gamma (\alpha) } \left[ \frac{x^{\alpha + 1}e^{-\beta x}}{- \beta } \right]^{\infty}_0 $$

$$ – \frac{\beta^{\alpha}}{\Gamma (\alpha) } \int_{0}^{\infty} \frac{(\alpha + 1)}{- \beta} x^{\alpha}e^{-\beta x} dx $$

第1項目は0になるため、
$$ = \frac{ (\alpha + 1) \beta^{\alpha-1}}{\Gamma (\alpha) } \int_{0}^{\infty} x^{\alpha}e^{-\beta x} dx $$
となり、ガンマ関数の性質より、\( \Gamma (\alpha + 1) = \alpha \Gamma (\alpha) \)から、
$$ = \alpha \beta^{-2} (\alpha + 1) \int_{0}^{\infty} \frac{\beta^{\alpha + 1}}{\Gamma (\alpha + 1)}x^{\alpha} e^{-\beta x} dx $$
ガンマ分布の定積分が1であることから、
$$ = \alpha \beta^{-2} (\alpha + 1) $$
となる。

$$ V(X) = E(X^2) – \left[ E(X) \right]^2 $$
$$ = \alpha \beta^{-2} (\alpha + 1) – \frac{\alpha^{2}}{\beta^{2}} $$
$$ = \frac{\alpha}{\beta^2} $$

5.指数分布

ガンマ分布の特別なケースとして、\( \alpha = 1, \beta = \frac{1}{\theta} \) とした、
$$ f(x; \theta) = \frac{1}{\theta} e^{-\frac{1}{\theta}x} , x > 0, \theta > 0 $$
の形をした分布を指数分布と呼ぶ。これはポアソン過程における\( n=1 \)、つまり初めて事象が起こるまでの時間の分布を表している。加えて、無記憶性があるため、幾何分布の連続型版と考えることもできる。

無記憶性について

$$ P(X > x + h | X > x) = \frac{\int^{\infty}_{x+h} \frac{1}{\theta} e^{-\frac{y}{\theta}} dy }{ \int^{\infty}_{x} \frac{1}{\theta} e^{-\frac{y}{\theta}} dy } $$

$$ = \frac{ \left[ \left( – \frac{1}{\theta} \right) e^{-\frac{y}{\theta}} \right]^{\infty}_{x+h} }{ \left[ \left( – \frac{1}{\theta} \right) e^{-\frac{y}{\theta}} \right]^{\infty}_{x} } $$
$$ = \frac{-e^{-\frac{x+h}{\theta}}}{-e^{-\frac{x}{\theta}}} = e^{-\frac{h}{\theta}} $$
よって、確率密度がこれまでの時間には依存しないことがわかる。

ハザードレートについて

あるt時点での確率密度を、t時点までの確率分布を1から引いたもので割ったものをハザードレートとする。(事象が起きるまでの確率1単位あたりの確率密度)
$$ \frac{f(t; \theta)}{1 – F(t;\theta)} = \frac{ \frac{1}{\theta} e^{-\frac{1}{\theta} t} }{ \int^{\infty}_{t} \frac{1}{\theta} e^{-\frac{1}{\theta} x} dx } $$
$$ = \frac{ \frac{1}{\theta} e^{-\frac{1}{\theta} t} }{ \left[ -e^{\frac{1}{\theta}x} \right]^{\infty}_{t} } = \frac{1}{\theta} $$
よって、tによらず一定であることがわかる。これは製品の寿命の分布と考えるに際して、時間を通じて寿命をむかえる確率が一定という指数分布の性質はクレイジーなものと思われますね。

ここで、指数分布の平均や分散を計算する上で必要な極限について取り上げておきます。

\( xe^{-x} \)の極限について

\( x > 0 \)として、
$$ x < 2^{x} $$
とし、両辺に\( e^{-x} \)をかけると、
$$ x e^{-x} < 2^{x} e^{-x} $$

ここで右辺に関して極限を取ると、
$$ \lim_{x \to \infty} 2^{x} e^{-x} = \lim_{x \to \infty} \left( \frac{2}{e} \right)^x = 0 $$
となる。
よって、
$$ \lim_{x \to \infty} x e^{-x} = 0 $$
となる。

\( x^n e^{-x} \)の極限について

\( x > 0 \)として、
$$ f(x) = x^{n+1}e^{-x} $$
を考える。
これを微分すると、
$$ f'(x) = (n+1)x^{n}e^{-x} – x^{n+1}e^{-x} $$
$$ = x^ne^{-x}(n+1-x) $$
ここで、\( f'(x) = 0 \)とすると、
$$ x = n+1 $$
となり、最大値は
$$ f(n+1) = (n+1)^{n+1} e^{-(n+1)} $$
となる。そこで以下の不等式から

$$ 0 < x^ne^{-x} = x^{n+1}e^{-x}x^{-1} \le (n+1)^{n+1}e^{-(n+1)}x^{-x} $$

最後の式について極限を取ると、
$$ \lim_{x \to \infty} (n+1)^{n+1}e^{-(n+1)}x^{-x} = 0 $$
となり、はさみうちされるため、
$$ \lim_{x \to \infty} x^n e^{-x} = 0$$
となる。

指数分布の平均

$$ E(X) = \int_{0}^{\infty} x \frac{1}{\theta} e^{ – \frac{1}{ \theta } x } dx $$
部分積分の公式より、
$$ = \left[ x \frac{-\theta}{\theta} e^{- \frac{1}{\theta}x} \right]^{\infty}_{0} -\int_{0}^{\infty} \frac{1}{\theta} (-\theta) e^{-\frac{1}{\theta}x} dx $$
\( xe^{-x} \)の極限より、第1項目は0になるため、
$$ = \int_{0}^{\infty} e^{-\frac{1}{\theta}x}dx $$
また指数分布の定積分が1であることから、
$$ = \theta \int_{0}^{\infty} \frac{1}{\theta} e^{-\frac{1}{\theta}x}dx = \theta $$
となる。

指数分布の分散

$$ E(X^2) = \int_{0}^{\infty} x^2 \frac{1}{\theta} e^{ – \frac{1}{ \theta } x } dx $$
部分積分の公式より、
$$ = \left[ x^2 \frac{-\theta}{\theta} e^{- \frac{1}{\theta}x} \right]^{\infty}_{0} – \int_{0}^{\infty} \frac{1}{\theta} (-\theta) 2x e^{-\frac{1}{\theta}x} dx $$

\( x^n e^{-x} \)の極限より、第1項目は0になるため、
$$ = \int_{0}^{\infty} 2x e^{-\frac{1}{\theta}x} dx $$
再び、部分積分の公式より、
$$ = 2 \left[ x( – \theta )e^{ – \frac{1}{ \theta } x } \right]^{\infty}_{0} – 2 \int_{0}^{\infty} ( – \theta )e^{-\frac{1}{\theta}x}dx $$

\( x e^{-x} \)の極限より、第1項目は0になるため、
$$ = 2 \int_{0}^{\infty} (\theta)e^{-\frac{1}{\theta}x}dx $$
$$ = 2 \theta^2 \int_{0}^{\infty} \frac{1}{\theta} e^{-\frac{1}{\theta}x}dx $$
指数分布の定積分が1であることから、
$$ = 2\theta^2 $$

$$ V(X) = E(X^2) – \left[ E(X) \right]^2 $$
$$ = 2\theta^2 – \theta^2 = \theta^2 $$

6.ワイブル分布

$$ f(x;\alpha, \beta) = \beta \alpha x^{\alpha-1}e^{-\beta x^{\alpha}}, x >0, \alpha > 0, \beta > 0 $$

の確率密度に従う分布をワイブル分布と呼び、\( \alpha = 1 \)のときは指数分布になる。

ハザードレートについて

あるt時点での確率密度を、t時点までの確率分布を1から引いたもので割ったものをハザードレートとする。(事象が起きるまでの確率1単位あたりの確率密度)

$$ \frac{f(t; \alpha, \beta)}{1 – F(t;\alpha, \beta)} = \frac{\beta \alpha t^{\alpha – 1} e^{-\beta t^{\alpha}}}{\int_{t}^{\infty}\beta \alpha x^{\alpha – 1}e^{-\beta x^{\alpha}}dx} $$

ここで、分母についてのみ注目して、\( z = x^{\alpha} \)と変数変換すると、\( \frac{dx}{dz} = \frac{1}{\alpha} z^{\frac{1}{\alpha}-1} \)から、

$$ \int_{t}^{\infty}\beta \alpha x^{\alpha – 1}e^{-\beta x^{\alpha}}dx = \beta \int_{t^{\alpha}}^{\infty} e^{-\beta z}dz $$

$$ = \beta \left[ – \frac{1}{\beta} e^{-\beta z} \right]^{\infty}_{t^{\alpha}} = e^{\beta t^{\alpha}} $$
よってハザードレートは、
$$ \frac{f(t; \alpha, \beta)}{1 – F(t;\alpha, \beta)} = \beta \alpha t^{\alpha -1} $$
となる。これは時間とともにハザードレートが高くなるという性質を示しており、全ての対象に対して、寿命に関して時間に対して単調増加を想定するのは誤りだろう。

ワイブル分布の平均

$$ E(X) = \int_{-\infty}^{\infty} x \beta \alpha x^{\alpha-1}e^{-\beta x^{\alpha}} dx $$
\( z = \beta x^{\alpha} \)として変数変換すると、\( \frac{dx}{dz} = \frac{1}{\alpha} z^{\frac{1}{\alpha}-1} \beta ^{-\frac{1}{\alpha}} \)から、
$$ = \int_{-\infty}^{\infty} \left( \frac{z}{\beta} \right)^{\frac{1}{\alpha}} \beta \alpha \left( \frac{z}{\beta} \right)^{\frac{\alpha – 1}{\alpha}} e^{-z} \frac{1}{\alpha}z^{\frac{1}{\alpha}-1} \beta ^{-\frac{1}{\alpha}} dz $$
$$ = \beta^{-\frac{1}{\alpha}} \int_{-\infty}^{\infty} z^{\frac{1}{\alpha}} e^{-z}dz $$
ガンマ関数の定義より、
$$ = \beta^{-\frac{1}{\alpha}} \Gamma \left( \frac{1}{\alpha} + 1 \right) $$

ワイブル分布の分散

$$ E(X^2) = \int_{-\infty}^{\infty} x^2 \beta \alpha x^{\alpha-1}e^{-\beta x^{\alpha}} dx $$
\( z = \beta x^{\alpha} \)として変数変換すると、
$$ = \int_{-\infty}^{\infty} \left( \frac{z}{\beta} \right)^{\frac{2}{\alpha}} \beta \alpha \left( \frac{z}{\beta} \right)^{\frac{\alpha – 1}{\alpha}} e^{-z} \frac{1}{\alpha}z^{\frac{1}{\alpha}-1} \beta ^{-\frac{1}{\alpha}} dz $$
$$ = \beta^{-\frac{2}{\alpha}} \int_{-\infty}^{\infty} z^{\frac{2}{\alpha}} e^{-z}dz $$
ガンマ関数の定義より、
$$ = \beta^{-\frac{2}{\alpha}} \Gamma \left( \frac{2}{\alpha} + 1 \right) $$
となる。

$$ V(X) = E(X^2) – \left[ E(X) \right]^2 $$
$$ = \beta^{-\frac{2}{\alpha}} \Gamma \left( \frac{2}{\alpha} + 1 \right) – \left[ \beta^{-\frac{1}{\alpha}} \Gamma \left( \frac{1}{\alpha} + 1 \right) \right]^2 $$
$$ = \beta^{-\frac{2}{\alpha}} \left[ \Gamma \left( \frac{2}{\alpha} + 1 \right) – \Gamma^2 \left( \frac{1}{\alpha} + 1 \right) \right] $$

7.ベータ分布

$$ f(x;\alpha, \beta) = \frac{1}{B(\alpha, \beta)}x^{\alpha – 1}(1-x)^{\beta – 1} , 0 < x < 1, \alpha > 0, \beta > 0 $$

ただし、\( B(\alpha, \beta) = \int_{0}^{1} x^{\alpha – 1}(1-x)^{\beta – 1} dx \)とする。この確率密度に従う分布をベータ分布と呼び、(0, 1)上の確率現象などをモデリングする際に用いられる。

ベータ関数の復習

ベータ関数は任意の2つの正定数x,yに対して定義される2変数関数で、
$$ B(x, y) = \int_{0}^{1} t^{x – 1}(1-t)^{y – 1} dt $$
で表され、収束する。
証明は区間を\( (0,\frac{1}{2}] \)と\( [ \frac{1}{2}, 1) \)に分けた積分について、各項が収束することを示すことによりなされる。

・ベータ関数の非負性、対称性

$$ B(x,y) > 0 $$

$$ B(x,y) = B(y,x) $$

対称性の証明は\( s = 1 – t \)と変数変換することで容易に示される。

・ベータ関数を三角関数の積分に変数変換

$$ B(x, y) = 2 \int_{0}^{\frac{\pi}{2}} \sin ^{2x-1} \theta \cos ^{2y – 1} \theta d \theta $$

証明は、\( t = \sin^2 \theta \)として変数変換をすることで容易に示される。

・ベータ関数とガンマ関数の関係について

$$ B(x, y) = \frac{ \Gamma (x) \Gamma (y) }{ \Gamma (x+y) } $$

証明はガンマ関数の定義式について、\( t = s^2 \)と変数変換を行い、極座標の変換公式を用いてから、三角関数で表現できるベータ関数を用いることで示される。

ベータ分布の平均

$$ E(X) = \int_{0}^{1} x \frac{1}{B(\alpha, \beta)}x^{\alpha – 1}(1-x)^{\beta – 1} dx $$

ベータ関数とガンマ関数の関係、\( B(x, y) = \frac{ \Gamma (x) \Gamma (y) }{ \Gamma (x+y) } \)から、
$$ = \int_{0}^{1} \frac{ \Gamma (\alpha + \beta ) }{ \Gamma (\alpha) \Gamma (\beta) } x^{\alpha}(1-x)^{\beta – 1} dx $$
ガンマ関数の性質\( \Gamma ( x + 1 ) = x \Gamma (x) \)より、

$$ = \int_{0}^{1} \frac{ \Gamma (\alpha + 1 + \beta ) }{ \alpha + \beta } \frac{\alpha}{\Gamma (\alpha + 1) \Gamma (\beta)} x^{\alpha}(1-x)^{\beta – 1} dx $$

$$ = \frac{\alpha}{\alpha + \beta} \int_{0}^{1} \frac{ \Gamma (\alpha + 1 + \beta ) }{ \Gamma (\alpha + 1) \Gamma (\beta) } x^{\alpha}(1-x)^{\beta – 1} dx $$
ベータ分布の定義より、定積分が1になることから、
$$ = \frac{\alpha}{\alpha + \beta} $$

ベータ分布の分散

$$ E(X^2) = \int_{0}^{1} x^2 \frac{1}{B(\alpha, \beta)}x^{\alpha – 1}(1-x)^{\beta – 1} dx $$
$$ = \int_{0}^{1} \frac{ \Gamma (\alpha + \beta ) }{ \Gamma (\alpha) \Gamma (\beta) } x^{\alpha+1}(1-x)^{\beta – 1} dx $$
ガンマ関数の性質より、

$$ = \int_{0}^{1} \frac{ \Gamma (\alpha + 1 + \beta ) }{ \alpha + \beta } \frac{\alpha}{\Gamma (\alpha + 1) \Gamma (\beta)} x^{\alpha+1}(1-x)^{\beta – 1} dx $$

さらに、ガンマ関数の性質より、

$$ = \int_{0}^{1} \frac{ \Gamma (\alpha + 2 + \beta ) }{ (\alpha + \beta) (\alpha + \beta + 1) } \frac{\alpha(\alpha + 1)}{\Gamma (\alpha + 2) \Gamma (\beta)} x^{\alpha+1}(1-x)^{\beta – 1} dx $$

整理すると、

$$ = \frac{\alpha(\alpha + 1)}{(\alpha + \beta) (\alpha + \beta + 1)} \int_{0}^{1} \frac{ \Gamma (\alpha + 2 + \beta ) }{\Gamma (\alpha + 2) \Gamma (\beta)} x^{\alpha+1}(1-x)^{\beta – 1} dx $$

となる。ベータ分布の定義より、定積分が1になることから、
$$ = \frac{\alpha(\alpha + 1)}{(\alpha + \beta) (\alpha + \beta + 1)} $$
となる。

$$ V(X) = E(X^2) – \left[ E(X) \right]^2 $$
$$ = \frac{\alpha(\alpha + 1)}{(\alpha + \beta) (\alpha + \beta + 1)} – \left[ \frac{\alpha}{\alpha + \beta} \right]^2 $$
$$ = \frac{\alpha(\alpha + 1)}{(\alpha + \beta) (\alpha + \beta + 1)} – \frac{\alpha^2}{(\alpha + \beta)^2} $$
$$ = \frac{\alpha(\alpha+1)(\alpha+\beta) – \alpha^2(\alpha + \beta +1)}{(\alpha + \beta)^2(\alpha + \beta +1)} $$
$$ = \frac{\alpha \beta}{(\alpha + \beta)^2(\alpha + \beta +1)} $$

8.コーシー分布

$$ f(x;\theta) = \frac{1}{\pi} \frac{1}{1+(x-\theta)^2} , – \infty < x < \infty, – \infty < \theta < \infty $$

という確率密度に従う分布をコーシー分布と呼ぶ。裾が長く、大きい値や小さい値を取る確率がなかなか0に近づかない分布とされる。
また、コーシー分布の特徴としては平均や分散を持たないことがあげられる。身近な活用例としては、株価の分析に使われたり、ベイズ統計学における無情報事前分布として半コーシー分布というものが用いられることがある。

ここで、\( x – \theta = z \)として、コーシー分布の平均を考える。
$$ E(X) = \int_{-\infty}^{\infty} z \frac{1}{\pi} \frac{1}{1+z^2} dz $$
上限をb、下限をaとすると、
$$ = \frac{1}{\pi} \int_{a}^{b} \frac{z}{1+z^2} dz$$
$$ = \frac{1}{\pi} \left[ \frac{1}{2} \log (1+z^2) \right]^b_a $$
$$ = \frac{1}{\pi 2} \log (1+b^2) – \frac{1}{\pi 2} \log (1+a^2) $$
ここで、aとbについて極限をとると、

$$ \lim_{b \to \infty} \lim_{a \to \infty} \left( \frac{1}{\pi 2} \log (1+b^2) – \frac{1}{\pi 2} \log (1+a^2) \right) $$

となり、これは不定となり極限値は存在しない。よって、平均値は存在しないため、分散も同様に存在しない。各々の項が無限にランダムに向かうため特定の値に収束することはないというイメージをするとわかりやすい。

9.多項分布

$$P(X_1=x_1, X_2=x_2, \dots, X_k=x_k) \\
:= \frac{n!}{x_1! x_2! \dots x_k!} \theta_1^{x_1} \theta_2^{x_2} \dots \theta_k^{x_k} \\
x_1, x_2, \dots x_k = 0, 1, \dots , n \\
x_1 + x_2 + \dots + x_k = n
$$

\( X_i \)は事象\( A_i \)の起こった回数を表し、\( x_i \)はその実現値を表す。事象\( A_i \)は標本空間において試行した結果として生じる。その試行は独立でn回繰り返しなされる。\( \theta_i = P( A_i ) \)である。
多項分布の期待値と分散の計算のために周辺結合分布を用いる必要があるため、ここでは多項分布の周辺化を行う。

\( X_i \)と\( X_j \)の周辺化

$$ P(X_i=x_i, X_j=x_j) \\
= \sum _{ \underset{ < i,j > } { x_1,\dots , x_k } } \frac{n!}{x_1! \dots x_k!} \theta_1^{x_1} \theta_2^{x_2} \dots \theta_k^{x_k} $$

\( x_i \)と\( x_j \)の項はここでは合計の対象として除かれるので、以下のように書き表すことができる。ただし後半の表現が一見すると何も変化がないように見えるため注意されたい。

$$ = \frac{ \theta_i^{x_i} \theta_j^{x_j} }{ x_i!x_j! } \sum _{ \underset{ < i,j > } { x_1,\dots , x_k } } \frac{n!}{x_1! \dots x_k!} \theta_1^{x_1} \theta_2^{x_2} \dots \theta_k^{x_k} $$

$$ = \frac{ \theta_i^{x_i} \theta_j^{x_j} (1 – \theta_i – \theta_j)^{n- x_i – x_j} n! }{ (n – x_i – x_j)! x_i! x_j! } \sum _{ \underset{ < i,j > } { x_1,\dots , x_k } } \frac{(n – x_i – x_j )!}{x_1! \dots x_k!} \left ( \frac{\theta_1}{1-\theta_i-\theta_j} \right )^{x_1} \left ( \frac{\theta_2}{1-\theta_i-\theta_j} \right )^{x_2} \dots \left ( \frac{\theta_k}{1-\theta_i-\theta_j} \right )^{x_k} $$

後半の部分はiとjを除いた多項分布のため、1となることから、
$$ = \frac{ n!}{ (n – x_i – x_j)! x_i! x_j! } \theta_i^{x_i} \theta_j^{x_j} (1 – \theta_i – \theta_j)^{n- x_i – x_j} $$

これは3項分布に他ならない。つまり、多項分布の2つの変数の周辺分布は3項分布となる。

\( X_i \)の周辺化

先程と同様にして周辺化を行う。
$$ P(X_i=x_i) = \sum _{ \underset{ < i > } { x_1,\dots , x_k } } \frac{n!}{x_1! \dots x_k!} \theta_1^{x_1} \theta_2^{x_2} \dots \theta_k^{x_k} $$

$$ = \frac{ \theta_i^{x_i} }{ x_i!} \sum _{ \underset{ < i > } { x_1,\dots , x_k } } \frac{n!}{x_1! \dots x_k!} \theta_1^{x_1} \theta_2^{x_2} \dots \theta_k^{x_k} $$

$$ = \frac{ \theta_i^{x_i} (1 – \theta_i )^{n- x_i } n! }{ (n – x_i)! x_i!} \sum _{ \underset{ < i > } { x_1,\dots , x_k } } \frac{(n – x_i )!}{x_1! \dots x_k!} \left ( \frac{\theta_1}{1-\theta_i} \right )^{x_1} \left ( \frac{\theta_2}{1-\theta_i} \right )^{x_2} \dots \left ( \frac{\theta_k}{1-\theta_i} \right )^{x_k} $$

$$ = \frac{ n! }{ (n – x_i)! x_i!} \theta_i^{x_i} (1 – \theta_i )^{n- x_i } \\
= {}_n C _{x_i} \theta_i^{x_i} (1 – \theta_i )^{n- x_i } $$

これは二項分布であることから、多項分布の一つの変数についての周辺分布は2項分布となる。

多項分布の平均

多項分布の期待値は二項分布のものと同じになる。

$$ E(X_i) = \sum _{ x_i \geq 0 , x_i \leq n } x_i \frac{ n! }{ (n – x_i)! x_i!} \theta_i^{x_i} (1 – \theta_i )^{n- x_i } $$

$$ = \sum _{ x_i \geq 0 , x_i \leq n } x_i \frac{ n! }{ (x_i – 1)! \\\{(n-1) – (x_i – 1)\\\}! } \theta_i^{x_i} (1 – \theta_i )^{n- x_i } $$
$$ = n \theta_i \sum _{ x_i \geq 1 , x_i \leq n } x_i \frac{ (n-1)! }{ (x_i – 1)! \\{(n-1) – (x_i – 1)\\}! } \theta_i^{x_i-1} (1 – \theta_i )^{(n-1)- (x_i-1) } $$

$$ = n \theta_i $$

多項分布の分散

多項分布の分散は二項分布のものと同じになる。

$$ V(X) = E(X^2) – \left [ E(X) \right ]^2 \\
= E(X(X-1)) + E(X) – \left [ E(X) \right ]^2
$$

$$ E(X_i(X_i-1)) \\
= \sum _{ x_i \geq 0 , x_i \leq n } x_i(x_i – 1) \frac{ n! }{ (n – x_i)! x_i!} \theta_i^{x_i} (1 – \theta_i )^{n- x_i } $$

$$ = \sum _{ x_i \geq 0 , x_i \leq n } \frac{ n! }{ \\{(n – 2) – (x_i -2 ) \\}! (x_i – 2)!} \theta_i^{x_i} (1 – \theta_i )^{n- x_i } $$
$$ = n(n-1) \theta_i^2 \\\
\sum _{ x_i \geq 2 , x_i \leq n } x_i \frac{ (n-2)! }{ (x_i – 2)! \\{(n-2) – (x_i – 2)\\}! } \theta_i^{x_i-2} (1 – \theta_i )^{(n-2)- (x_i-2) } $$

$$= n(n-1) \theta_i^2$$

$$ V(X) = E(X(X-1)) + E(X) – \left [ E(X) \right ]^2 \\
= n(n-1) \theta_i^2 + n \theta_i – ( n \theta_i )^2 \\
= n \theta_i (1 – \theta_i)
$$

多項分布の共分散

$$
E(X_i X_j) = \sum _{ x_i,x_j \geq 0 \\\ x_i + x_j \leq n } x_i x_j \frac{n!}{x_i! x_j!(n-x_i-x_j)!} \theta_i^{x_i}\theta_j^{x_j}(1-\theta_i – \theta_j)^{n-x_i-x_j}
$$
$$ = \sum _{ x_i,x_j \geq 0 \\\ x_i + x_j \leq n } \frac{n!}{(x_i-1)! (x_j-1)![(n-2)-(x_i-1)-(x_j-1)]!} \theta_i^{x_i}\theta_j^{x_j}(1-\theta_i – \theta_j)^{n-x_i-x_j} $$
$$ = n(n-1) \theta_i \theta_j \sum _{ x_i,x_j \geq 1 \\\ x_i + x_j \leq n } \frac{(n-2)!}{(x_i-1)! (x_j-1)![(n-2)-(x_i-1)-(x_j-1)]!} \theta_i^{x_i-1}\theta_j^{x_j-1}(1-\theta_i – \theta_j)^{(n-2)-(x_i-1)-(x_j-1)} $$

$$ = n(n-1)\theta_i \theta_j $$

$$ Cov(X_i, X_j) = E( X_i X_j ) – E(X_i) E(X_j) \\
= n(n-1)\theta_i \theta_j – n\theta_i n\theta_j \\
= – n \theta_i \theta_j
$$

多項分布の条件付き分布

\( X_j = x_j \)のときの\( X_i \)の条件付き分布は\( P(X_i = x_i | X_j = x_j) = \frac{P(X_i \cap X_j)}{P(X_j)} \)より、

$$ P(X_i = x_i | X_j = x_j) = \frac{ \frac{n!}{x_i!x_j!(n-x_i-x_j)!} \theta_i^{x_i} \theta_j^{x_j} (1-\theta_i – \theta_j)^{n-x_i-x_j} }{ \frac{n!}{x_j!(n-x_j)!} \theta_j^{x_j} (1- \theta_j)^{n-x_j} } $$
$$ = \frac{(n-x_j)!}{x_i!(n-x_i-x_j)!} \theta_i^{x_i} \left ( \frac{1}{1-\theta_j} \right )^{x_j – n} \left ( 1 – \theta_i – \theta_j \right )^{n-x_i-x_j} $$
$$ = \frac{(n-x_j)!}{x_i!(n-x_i-x_j)!} \left ( \frac{\theta_i}{1-\theta_j} \right )^{x_i} \left ( \frac{1}{1-\theta_j} \right )^{n- x_i -x_j} \left ( 1 – \theta_i – \theta_j \right )^{n-x_i-x_j} $$
$$ = \frac{(n-x_j)!}{x_i!(n-x_i-x_j)!} \left ( \frac{\theta_i}{1-\theta_j} \right )^{x_i} \left ( \frac{ 1 – \theta_i – \theta_j }{1-\theta_j} \right )^{n- x_i -x_j} $$
$$ = \frac{(n-x_j)!}{x_i!(n-x_i-x_j)!} \left ( \frac{\theta_i}{1-\theta_j} \right )^{x_i} \left ( 1 – \frac{\theta_i}{1-\theta_j} \right )^{n- x_i -x_j} $$

$$ = {}_{n-x_j} C _{x_i} \hat \theta_i ^{x_i}(1 – \hat \theta_i)^{(n-x_j)-x_i} $$

$$ \hat \theta_i = \frac{\theta_i}{1-\theta_j} $$

よってこれは二項分布であるから、多項分布の\( X_j = x_j \)のときの\( X_i \)の条件付き分布は二項分布に従うことがわかる。

10.ディリクレ分布

ディリクレ分布は\( \sum_{i=1}^{m} \theta_i = 1, \theta_i \geq 0 \)を満たす\( \theta_1, \dots , \theta_m \)の集合からなる。
ディリクレ分布の確率密度関数は以下の形をとる。

$$ f(\theta_1, \dots , \theta_m) = \frac{\Gamma (A)}{ \prod_{i=1}^{m} \Gamma (\alpha_i) } \prod_{i=1}^{m} \theta_i^{\alpha_i – 1} \\
A = \sum_{i=1}^{m} \alpha_i \\
\alpha_i > 0, \quad i = 1, \dots , m
$$

\( m = 2 \)のとき、ディリクレ分布はベータ分布の特殊型になる。

$$ f(\theta_1, \theta_2) = \frac{ \Gamma ( \alpha_1 + \alpha_2 ) }{ \Gamma (\alpha_1) \Gamma (\alpha_2) } \theta_1^{\alpha_1 – 1} \theta_2^{\alpha_2 – 1} $$

ディリクレ分布の平均

$$ E(\theta_j) = \int \dots \int \theta_j \frac{\Gamma (A)}{ \prod_{i=1}^{m} \Gamma (\alpha_i) } \prod_{i=1}^{m} \theta_i^{\alpha_i – 1} d\theta_1 \dots d\theta_{m-1} $$
$$ = \frac{\Gamma(A)}{\Gamma(A+1)} \frac{\Gamma(\alpha_j + 1)}{\Gamma(\alpha_j)} \int \dots \int \frac{\Gamma (A+1)}{ \prod_{i=1}^{m} \Gamma (\alpha_i^{\prime}) } \prod_{i=1}^{m} \theta_i^{\alpha_i^{\prime} – 1} d\theta_1 \dots d\theta_{m-1} $$

ここでは\(
\begin{cases}
\alpha_i^{\prime} = \alpha_i \quad i \neq j \\
\alpha_j^{\prime} = \alpha_j + 1 \quad i = j \\
\end{cases}
\)となる。

積分以降はディリクレ分布の定義より1なので、以下のように表される。

$$ = \frac{\Gamma(A)}{\Gamma(A+1)} \frac{\Gamma(\alpha_j + 1)}{\Gamma(\alpha_j)} $$
ガンマ関数の定義より、

$$ = \frac{\Gamma(A)}{\Gamma(A)A } \frac{\Gamma(\alpha_j)\alpha_j}{\Gamma(\alpha_j)} $$
となることから、
$$ = \frac{\alpha_j}{A} $$
となる。

ディリクレ分布の分散

分散は以下の式より求まる。
$$ V(\theta_j) = E \left( \theta_j^2 \right) – \left [ E(\theta_j) \right ]^2 $$

$$ E(\theta_j^2) = \int \dots \int \theta_j^2 \frac{\Gamma (A)}{ \prod_{i=1}^{m} \Gamma (\alpha_i) } \prod_{i=1}^{m} \theta_i^{\alpha_i – 1} d\theta_1 \dots d\theta_{m-1} $$
$$ = \frac{\Gamma(A)}{\Gamma(A+2)} \frac{\Gamma(\alpha_j + 2)}{\Gamma(\alpha_j)} \int \dots \int \frac{\Gamma (A+2)}{ \prod_{i=1}^{m} \Gamma (\alpha_i^{\prime}) } \prod_{i=1}^{m} \theta_i^{\alpha_i^{\prime} – 1} d\theta_1 \dots d\theta_{m-1} $$

$$ = \frac{\Gamma(A)}{\Gamma(A+1)(A+1)} \frac{\Gamma(\alpha_j + 1)(\alpha_j + 1)}{\Gamma(\alpha_j)} $$

$$ = \frac{\Gamma(A)}{\Gamma(A)A(A+1)} \frac{\Gamma(\alpha_j)\alpha_j(\alpha_j + 1)}{\Gamma(\alpha_j)} $$

$$ = \frac{\alpha_j (\alpha_j + 1)}{A(A+1)} $$

$$ V(\theta_j) = E \left( \theta_j^2 \right) – \left [ E(\theta_j) \right ]^2 \\
= \frac{\alpha_j (\alpha_j + 1)}{A(A+1)} – \left ( \frac{\alpha_j}{A} \right )^2 \\
= \frac{\alpha_j}{A(A+1)} – \frac{\alpha_j^2}{A^2(A+1)}
$$

多変量正規分布

追記予定

多変量正規分布の平均

追記予定

多変量正規分布の分散

追記予定

参考文献

[1] 鈴木・山田 (1996), 『数理統計学―基礎から学ぶデータ解析』, 内田老鶴圃
[2] 原隆 , “微分積分続論 SII-15 クラス” , 九州大学の講義資料
[3] 高校数学の美しい物語 , “ガウス積分の公式の2通りの証明”
[4] 高校数学の美しい物語 , “偶関数と奇関数の意味,性質などまとめ”
[5] 小杉考司 (2015) ,”Cauchy分布について(ベイズ塾例会資料)2015.07.26″, slideshare
[6] MAS3301 Bayesian Statistics

[数理統計学]離散型確率分布の期待値と分散の導出まとめ

はじめに

先日、知人が社労士の試験を受けていたのですが、アプリで問題を色々と解けるものがあるらしく、少しの課金で通勤電車などでの勉強が捗るそうです。統計検定に関してもそのようなものがあればいいなとは思いますが、そんなものを作る気力はないので、よくある分布とその平均と分散に関する導出についてひたすら記していきます。これで少なくとも自分の通勤時の学習が捗ると思われます。なによりも、歳をとっても最低限の数式を扱うスキルがあれば思い出せるレベルで残すことも意識したいですね。

※PC版でないと数式が切れてしまうので、SP版での閲覧はおすすめしません。

今回登場する離散分布

『数理統計学―基礎から学ぶデータ解析』という本に登場するものを式を補いながら紹介します。

  • 1.一様分布
  • 2.ベルヌーイ分布
  • 3.二項分布
  • 4.ポアソン分布
  • 5.超幾何分布
  • 6.幾何分布
  • 7.負の二項分布

1.一様分布

確率変数が有限個の値をとり、それぞれの確率が$$\frac{1}{n}$$である分布。

一様分布の平均

$$E(X)=\frac{1}{n}\sum_{i=1}^{n}X_{i}$$

一様分布の分散

$$V(X)=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\bar X)^2$$

2.ベルヌーイ分布

等確率の原理に従うとして、あるトラックの中にある「ちくわぶ」を取り出していった場合、確率θで「ちくわぶ」は穴が空いておらず、確率(1-θ)で穴が空いているとする。穴の空いていない「ちくわぶ」の割合はベルヌーイ分布に従う。

ベルヌーイ分布の平均

$$E(X)=1 \times \theta + 0 \times (1-\theta ) \ $$
$$ = \theta $$

ベルヌーイ分布の分散

$$V(X)=(1-\theta)^2 \times \theta + (0 – \theta)^2 \times (1-\theta ) \\ $$
$$ = \left[ (1-\theta) \times \theta + \theta^2 \right] \times (1-\theta) $$

$$ = \theta(1-\theta) $$

3.二項分布

一つの事象(「ちくわぶ」に穴が空いているかどうか)に関する、確率θからなるベルヌーイ分布に従う確率変数列を
$$X_1 , X_2, \dots ,X_n$$
として、それらの和(n回試行した際の「ちくわぶ」に穴が空いていない個数)を
$$X = X_1 + X_2 + \dots + X_n$$
とする。その和の従う確率分布は二項分布と呼ばれ、
$$P(X=x):= {}_n C_x \theta^x(1-\theta)^{n-x}$$
と表される。

二項分布の平均

$$E(X)=\sum_{x=0}^{n}x \times {}_n C_x \theta^x(1-\theta)^{n-x} $$

x=0のときは0なので、

$$ = \sum_{x=1}^{n} \frac{n!}{(n-x)!(x-1)!}\theta^x(1-\theta)^{n-x}$$
$$ = n \theta \sum_{x=1}^{n} \frac{(n-1)!}{(n-x)!(x-1)!}\theta^{x-1}(1-\theta)^{n-x}$$
$$ = n \theta \sum_{x=1}^{n} \frac{(n-1)!}{\left[ (n-1) – (x-1) \right]!(x-1)!}$$
$$ \times \theta^{x-1}(1-\theta)^{ (n-1) – (x-1) } $$
ここで
$$y = x-1$$
とおくと、

$$ = n \theta \sum_{y=0}^{n-1} {}_{n-1} C_y \theta^y(1-\theta)^{n-1-y}$$
二項分布の和は1であるという性質を使えば、
$$ = n \theta $$

二項分布の分散

$$ E [ X(X-1) ] = \sum_{x=0}^{n} x(x-1) $$
$$ \times {}_n C_x \theta^x(1-\theta)^{n-x} $$

x = 1, 2のときは0なので、

$$= \sum_{x=2}^{n} \frac{n!}{(n-x)!(x-2)!} \theta^x(1-\theta)^{n-x} $$
$$= \sum_{x=2}^{n} \frac{n!}{\left[ (n-2) – (x-2) \right]!(x-2)!} \theta^x(1-\theta)^{(n-2) – (x-2) } $$
$$= n(n-1)\theta^2 \sum_{x=2}^{n} \frac{n!}{\left[ (n-2) – (x-2) \right]!(x-2)!} \theta^{x-2}(1-\theta)^{(n-2) – (x-2) } $$

ここで
$$y = x -2$$
とおくと、

$$= n(n-1)\theta^2 \sum_{y=0}^{n-2} {}_{n-2} C_y \theta^y(1-\theta)^{(n-2) – y } $$

となり、二項分布の和が1であることから、
$$= n(n-1)\theta^2$$
となる。

$$V(X)=E(X^2)-[ E(X) ]^2$$
$$ =E[X(X-1)] + E(X) -[ E(X) ]^2$$
$$ =n(n-1)\theta^2 + n \theta -[ n \theta ]^2$$
$$ =n\theta(1-\theta)$$

また、この式より、二項分布に関してはθの絶対値が1よりも小さいことから、平均よりも分散の方が小さくなることがわかる。

4.ポアソン分布

二項分布において、平均値がλとおく。つまり、
$$n\theta = \lambda$$
λが一定でnが非常に大きいケースにおいて、その確率分布はポアソン分布に従う。λが一定でnが非常に大きいケースというのは、θが非常に小さいということになり、その事象が滅多に起きないような対象に対して扱われることが多い。

ポアソン分布の導出は二項分布からネイピア数を用いて導出する方法や、微分方程式を用いた方法などがある。

ポアソン分布の導出(二項分布からver)

$${}_{n} C_x \theta^x(1-\theta)^{n-x}$$
$$ = \frac{n(n-1) \dots (n-x+1)}{x!} \theta^x(1-\theta)^{n-x}$$

ここで、
$$\theta = \frac{\lambda}{n}$$
とすると

$$ = \frac{1\left(1-\frac{1}{n}\right) \dots \left(1-\frac{x-1}{n}\right)}{x!} \left(\frac{\lambda}{n}\right)^x\left(1-\frac{\lambda}{n}\right)^{n-x}n^x$$
$$ = \frac{1\left(1-\frac{1}{n}\right) \dots \left(1-\frac{x-1}{n}\right)}{x!} \lambda^x\left(\left(1-\frac{\lambda}{n}\right)^{-\frac{n}{\lambda}}\right)^{-\lambda} \left(1-\frac{\lambda}{n}\right)^{-x}$$

nを無限大にすると、分母にのみnがある項は0になり、また、ネイピア数の定義から
$$ \lim_{n\to\infty}{}_{n} C_x \theta^x(1-\theta)^{n-x} = \frac{\lambda^x e^{-\lambda}}{x!}$$
となる。
$$P(X=x):= \frac{\lambda^x e^{-\lambda}}{x!}$$
この確率分布がポアソン分布と呼ばれる。(x=0,1,2,…)

実際に、確率分布の性質を確かめると、

$$ \sum_{x=0}^{\infty} P(X=x) = \sum_{x=0}^{\infty} \frac{\lambda^x e^{-\lambda}}{x!} $$
$$ = e^{-\lambda} \sum_{x=0}^{\infty} \frac{\lambda^x}{x!} $$
指数の和の公式より、
$$ = e^{-\lambda} e^{\lambda} = 1 $$
よって、確率分布であることがわかる。

ポアソン分布の導出(微分方程式ver)

abrahamcowさんのブログにまさにそれが書かれていました。図も載っていて詳しくて良いです。
他にも、英語の文献だとA Quick Way to See that the Poisson Distribution is the Appropriate Mathematical Formulation for a Counting Process with Constant Rate and Intensityなどがあります。

これらの文献を見ていただければ導出について理解できると思います。

以下、これらの文献を参考にして式展開したものを記しておきます。

まず、3つの仮定をおきます。

  • 1.互いに背反する期間では現象が起こる確率は独立。
    つまり、X(t)とX(t+h)-X(t)は独立。
  • 2.非常に小さい区間で現象が起こる確率はその区間の長さに比例する。
    $$P\{X(t+h) – X(t) = 1 \} = \lambda h + o(h)$$
    ここでo(h)は
    $$\lim_{n\to 0} \frac{o(h)}{h}=0$$
    となる関数とする。また、λ > 0とする。
  • 3.非常に小さい区間で現象が2回以上起こる確率はその区間で現象が1回起こる確率に比べて、無視できるほど小さい。
    $$P\{X(t+h) – X(t) = k\} = o(h) $$
    ただし、k > 1とする。

ここで、X(t)がrになる確率を以下の表現で表す。
$$P_r (t) = P\{ X(t) = r \} $$

事象
$$X(t+h)=r$$
は以下の3つに場合分けできる。

  • 1.$$ \{ X(t) = r,X(t+h) – X(t) = 0 \}$$
    この
    $$X(t+h) – X(t) = 0$$
    は仮定2と仮定3の式を1から差っ引いたものによって確率を計算できる。また、仮定1から互いに背反する期間では現象が起こる確率は独立のため、同時確率は各々の積で表される。

$$P\{ X(t) = r \} (1 – \lambda h – o(h) – o(h) )$$

  • 2.$$ \{ X(t) = r-1,X(t+h) – X(t) = 1 \}$$
    仮定2より、
    $$P\{ X(t) = r-1 \} (\lambda h + o(h))$$
    によって表される。

  • 3.

$$ \\{ X(t+h) = r-k,X(t+h) – X(t) = k \\} , k > 1 $$

仮定3より、
$$P\{ X(t+h) = r-k \} o(h)$$
によって表される。

これらは期間が違うため、仮定より互いに背反であることから、3つに場合分けの確率の和で表現する。
$$P_r(t+h) = P\{ X(t) = r \} $$
$$= P\{ X(t) = r \} (1 – \lambda h – o(h) – o(h) ) $$

$$+ P\{ X(t) = r-1 \} (\lambda h + o(h)) + P\{ X(t+h) = r-k \} o(h)$$

$$ = P_r(t)(1-\lambda h) + P_{r-1}(t)\lambda h + o(h) $$

$$P_r(t+h) – P_r(t) = -\lambda h P_r(t) + P_{r-1}(t)\lambda h + o(h) $$
$$\frac{P_r(t+h) – P_r(t)}{h} = -\lambda P_r(t) + \lambda P_{r-1}(t) + \frac{o(h)}{h} $$

hの極限をとると、

$$\lim_{h\to 0} \frac{P_r(t+h) – P_r(t)}{h} = -\lambda P_r(t) + \lambda P_{r-1}(t) $$

となる。

これらは
r = 0のとき
$$\frac{P_0(t)}{dt} = – \lambda P_0(t)$$

r ≠ 0のとき

$$\frac{P_r(t)}{dt} = -\lambda P_r(t) + \lambda P_{t-1}(t), r = 1,2,\dots$$

となる微分方程式として表される。

まず、r=0の場合の微分方程式を解くと、

$$\frac{P_0(t)}{dt} = – \lambda P_0(t) \Rightarrow \frac{d P_0(t)}{P_0(t)} = -\lambda dt $$
両辺積分すると、
$$ \int \frac{d P_0(t)}{P_0(t)} = -\lambda \int dt $$

$$ \log P_0(t) = -\lambda t + C$$
(Cは積分定数)

$$P_0(t) = e^{-\lambda t +C}$$
Cは任意なので、
$$e^C=C$$
として、
$$P_0(t) = Ce^{-\lambda t}$$
となる。
t = 0で1件も発生しない確率は1なので、
$$P_0(0)=C=1$$
から、C = 1であることがわかり、
$$P_0(t) = e^{-\lambda t}$$
となる。

ここで、定数係数の1階線形微分方程式の解の公式は、
$$y^{\prime} + p(x)y = q(x)$$
に対して、
$$y = e^{-p(x)x}\left[ \int q(x) e^{p(x)x}dx + C \right]$$
で表される。
この解の公式を今回のr=1の場合の微分方程式に当てはめると、
$$y^{\prime} = P_1(t) , p(x) = -\lambda , q(x) = \lambda P_0(t) $$
であるから、
$$P_1(t) = e^{-\lambda t}\left[ \int \lambda P_0(t) e^{\lambda t} dt + C \right] $$
となる。ここに、
$$P_0(t)$$
を代入すると、
$$ = e^{-\lambda t}\left[ \int \lambda e^{-\lambda t} e^{\lambda t} dt \right] $$
$$ = e^{-\lambda t}\left[ \int \lambda dt + C \right] $$
$$ = e^{-\lambda t} \lambda t + e^{-\lambda t}C $$
となる。
t=0においてX(0)=1の確率は0なので、C=0となり、
$$P_1(t) = e^{-\lambda t} \lambda t $$

同様にしてr=2の場合、
$$P_2(t) = \frac{e^{-\lambda t \left( \lambda t \right)^2}}{2} $$
r=3の場合、
$$P_3(t) = \frac{e^{-\lambda t \left( \lambda t \right)^3}}{2 \times 3} $$
となる。

ここで、任意の整数r=kにおいて、
$$P_k(t) = \frac{e^{\lambda t} \left( \lambda t \right)^k }{k!}$$
が成り立つとする。

r=k+1のとき、
$$P_{k+1}(t) = \frac{e^{\lambda t} \left( \lambda t \right)^{k+1} }{(k+1)!}$$
である。
これはポアソン分布であり、定数係数の1階線形微分方程式よりポアソン分布が導出できることが示された。

ポアソン分布の平均

$$E(X)=\sum_{x=0}^{\infty}x \frac{e^{-\lambda}\lambda^x}{x!}$$
$$ =\sum_{x=0}^{\infty} \frac{e^{-\lambda}\lambda^x}{(x-1)!}$$
$$ =\lambda \sum_{x=1}^{\infty} \frac{e^{-\lambda}\lambda^{x-1}}{(x-1)!} = \lambda $$

ポアソン分布の分散

$$E(X(X-1)) = \sum_{x=0}^{\infty}x(x-1) \frac{e^{-\lambda}\lambda^x}{x!}$$
$$ = \sum_{x=0}^{\infty} \frac{e^{-\lambda}\lambda^x}{(x-2)!}$$
$$ = \lambda^2 \sum_{x=2}^{\infty} \frac{e^{-\lambda}\lambda^{x-2}}{(x-2)!} =\lambda^2 $$

$$V(X)=E(X^2)-[ E(X) ]^2$$
$$ =E[X(X-1)] + E(X) -[ E(X) ]^2$$
$$ = \lambda^2 + \lambda – \lambda^2 = \lambda $$

ポアソン分布の再生性

XとYは独立にそれぞれポアソン分布λ1とλ2に従うとする。
$$P(Z=z) = \sum_{x=0}^{z} P(X=x, Y=z-x) $$
$$ = \sum_{x=0}^{z} P(X=x)P(Y=z-x) $$

$$ = \sum_{x=0}^{z} \frac{e^{-\lambda_1} \lambda_1^x }{x!} \frac{e^{-\lambda_2} \lambda_2^{z-x} }{(z-x)!} $$
$$ = \frac{e^{-(\lambda_1 + \lambda_2 )}(\lambda_1 + \lambda_2)^z}{z!} \sum_{x=0}^{z}\frac{z!}{x!(z-x)!}\left( \frac{\lambda_1}{\lambda_1+\lambda_2} \right)^x \left( \frac{\lambda_2}{\lambda_1+\lambda_2} \right)^{z-x} $$

二項分布の和は1になるので以下のように表される。
$$ = \frac{e^{-(\lambda_1 + \lambda_2 )}(\lambda_1 + \lambda_2)^z}{z!} $$
このように表されることを再生性があると呼ぶ。

5.超幾何分布

N個の製品からなる仕切りのなかで、M個の不良品があるとする。
非復元抽出でn個を抽出した際の、不良品の数をXとした際の確率分布は以下のように表される。

$$P(X=x)=\frac{{}_M \mathrm{C} _x \times {}_{N-M} \mathrm{C} _{n-x}}{{}_N \mathrm{C} _n} $$

xは非負の整数で

$$max \{ 0, n -(N-M) \} \leq x \leq min \{ n, M \}$$

に従うものとする。

ここで製品の数Nが標本nに比べて十分に大きい場合を考える。
$$\frac{{}_M \mathrm{C} _x \times {}_{N-M} \mathrm{C} _{n-x}}{{}_N \mathrm{C} _n} $$

$$ = \frac{M!}{x!\left( M-x \right)!} \times \frac{\left( N – M \right)! }{\left( n-x\right)! \left[ N-M -\left( n-x\right) \right]!} \times \frac{n! \left(N-n\right)!}{N!} $$
$$ = \frac{\left[ M(M-1)\dots (M-x+1)\right] }{x!(n-x)! \times \left[ N(N-1)\dots (N-n+1) \right]}$$
$$ \times \left[ (N-M)(N-M-1)\dots (N-M-n+x+1)\right]\times n! $$

組み合わせを使ってまとめると、分母と分子の数はn個であるから、分母分子にn個だけ1/Nを掛けると以下のようになる。
$$ = \frac{{}_n \mathrm{C} _x \times \frac{M}{N} \left( \frac{M}{N} – \frac{1}{N} \right) \dots \left( \frac{M}{N} – \frac{x-1}{N} \right) }{1 \left( 1 – \frac{1}{N} \right) \dots \left( 1 – \frac{n-1}{N} \right) } $$

$$ \times \left( \frac{N-M}{N} \right) \times \left( \frac{N-M}{N} – \frac{1}{N} \right) \dots \left( \frac{N-M}{N} – \frac{n-x-1}{N} \right) $$

ここで
$$ \theta = \frac{M}{N}$$
から、
$$ M = N\theta$$
として代入すると、

$$ = \frac{{}_n \mathrm{C} _x \times \theta \left( \theta – \frac{1}{N} \right) \dots \left( \theta – \frac{x-1}{N} \right) }{1 \left( 1 – \frac{1}{N} \right) \dots \left( 1 – \frac{n-1}{N} \right) } $$
$$ \times \left( 1-\theta \right) \times \left( (1-\theta) – \frac{1}{N} \right) \dots \left( (1-\theta) – \frac{n-x-1}{N} \right) $$

となる。
ここで、Nの極限を取ると以下のように、分母が全て1になり、以下のように表される。

$$ \lim_{N\to\infty}\frac{{}_M \mathrm{C} _x \times {}_{N-M} \mathrm{C} _{n-x}}{{}_N \mathrm{C} _n} $$
$$ = {}_n \mathrm{C} _x \theta^x \left( 1 – \theta \right)^{n-x} $$

これは二項分布と同じ形となる。
つまり、超幾何分布において、製品の数が十分に大きいと二項分布と同じになる。

超幾何分布の平均

$$ E(X) = \sum_{x=0}^{n} x \frac{ {}_M \mathrm{C} _x \times {}_{N-M} \mathrm{C} _{n-x}}{{}_N \mathrm{C} _n} $$

$$ = \frac{1}{ {}_N \mathrm{C} _n } \sum_{x=0}^{n} x \frac{M!}{(M-x)!x!} $$
$$ \times \frac{(N-M)!}{ \left[ N-M-(n-x) \right]! (n-x)!} $$

$$ = \frac{M}{ {}_N \mathrm{C} _n} \sum_{x=0}^{n} \frac{(M-1)!}{\left[ (M-1)-(x-1) \right]!(x-1)!} $$

$$ \times \frac{\left[ (N-1)-(M-1) \right]!}{\left[ (N-1)-(M-1)-\left[(n-1)-(x-1)\right] \right]! \left[ (n-1)-(x-1) \right]!}$$

$$ = \frac{M}{ {}_N \mathrm{C} _n} $$
$$ \times \sum_{x=0}^{n} {}_{M-1} \mathrm{C} _{x-1} \times {}_{(N-1) – (M-1)} \mathrm{C} _{(n-1)-(x-1)} $$

$$ = \frac{M}{ \frac{N!}{(N-n)!n!} } $$
$$ \times \sum_{x=0}^{n} {}_{M-1} \mathrm{C} _{x-1} \times {}_{(N-1) – (M-1)} \mathrm{C} _{(n-1)-(x-1)} $$

$$ = \frac{nM}{ \frac{(N-1)!N}{\left[ (N-1)-(n-1)\right]!(n-1)!} } $$
$$ \times \sum_{x=0}^{n} {}_{M-1} \mathrm{C} _{x-1} \times {}_{(N-1) – (M-1)} \mathrm{C} _{(n-1)-(x-1)}$$

$$ = \frac{nM}{N} $$
$$ \sum_{x=0}^{n} \frac{ {}_{M-1} \mathrm{C} _{x-1} \times {}_{(N-1) – (M-1)} \mathrm{C} _{(n-1)-(x-1)} }{ {}_{N-1} \mathrm{C} _{n-1} } $$

超幾何分布の和は確率分布のため、1になることから、
$$ = \frac{nM}{N}$$
となる。これは二項分布の平均と同じとなる。

いちいち計算するのが面倒なので、二項係数の公式をここに載せておきます。

$$ {}_n \mathrm{C} _r = {}_n \mathrm{C} _{n-r} $$
$$r {}_n \mathrm{C} _r = n {}_{n-1} \mathrm{C} _{r-1}$$
$$ {}_n \mathrm{C} _r = {}_{n-1} \mathrm{C} _{r} + {}_{n-1} \mathrm{C} _{r-1}$$

超幾何分布の分散

$$E(X(X-1))= \sum_{x=0}^{n} x(x-1) $$
$$ \times \frac{ {}_M \mathrm{C} _x \times {}_{N-M} \mathrm{C} _{n-x}}{{}_N \mathrm{C} _n} $$

二項係数の公式より、

$$ = \sum_{x=0}^{n} (x-1)$$
$$ \times \frac{ M \times {}_{M-1} \mathrm{C} _{x-1} \times {}_{N-M} \mathrm{C} _{n-x} \times n }{ N \times {}_{N-1} \mathrm{C} _{n-1} } $$

$$ = \sum_{x=0}^{n} M(M-1) \times n(n-1) $$
$$ \times \frac{ {}_{M-2} \mathrm{C} _{x-2} \times {}_{(N-2)-(M-2)} \mathrm{C} _{(n-2)-(x-2)} }{ N(N-1) \times {}_{N-2} \mathrm{C} _{n-2}} $$

超幾何分布の和が1であることから、
$$ n(n-1) \frac{M(M-1)}{N(N-1)} $$

$$V(X)=E(X^2)-[ E(X) ]^2$$
$$ =E[X(X-1)] + E(X) -[ E(X) ]^2$$
$$ = n(n-1) \frac{M(M-1)}{N(N-1)}+ \frac{nM}{N} – \frac{n^2M^2}{N^2} $$
$$ = \frac{nM}{N} \left[ (n-1) \frac{M-1}{N-1} + 1 – \frac{nM}{N} \right] $$
$$ = \frac{nM}{N} \left( 1 – \frac{M}{N} \right) \left( \frac{N-n}{N-1} \right) $$

超幾何分布の平均は二項分布と同じであったが、分散に関しては、
$$ \left( \frac{N-n}{N-1} \right) $$
だけ異なる。

6.幾何分布

復元抽出で、不良品が出るまで検査した際の観測された良品の数Xの確率分布は以下の幾何分布に従う。

$$P(X=x)=\theta ( 1-\theta)^x $$

幾何分布は無記憶性という性質を持っており、それは幾何分布のみが持つとされている。無記憶性はある時点から先に時点を進めた際に、過去の時点の影響が残らないことを指している。

$$P(X=x+h | X \geq x) = \frac{P(X=x+h)}{ \sum_{y=x}^{\infty} P(X=y) }$$

$$ = \frac{ \theta ( 1-\theta)^{x+h} }{ \sum_{y=x}^{\infty} \theta ( 1-\theta)^y } $$
$$ = \frac{ \theta ( 1-\theta)^{h} }{ \sum_{y=x+1}^{\infty} \theta ( 1-\theta)^y } $$
超幾何分布の和は1であるから、
$$ = \theta ( 1-\theta)^h $$
$$ = P(X=h) $$

幾何分布の平均

$$ E(X) = \sum_{x=0}^{\infty} x \theta ( 1-\theta)^x $$
$$ = \sum_{x=1}^{\infty} x \theta ( 1-\theta)^x $$
$$ = \theta \sum_{x=1}^{\infty} x ( 1-\theta)^x $$
$$ = \theta \frac{d \sum_{x=1}^{\infty} (1-\theta)^x}{d(1-\theta)}\times (1-\theta) $$
無限等比級数の和の公式より、
$$ = \theta \frac{d \left( \frac{1-\theta}{1-(1-\theta)} \right) }{d(1-\theta)}\times (1-\theta) $$
$$ = \theta \frac{d \frac{1-\theta}{\theta}}{d \theta} \times \frac{d \theta}{d(1-\theta)}\times (1-\theta)$$
$$ = \theta \frac{-\theta – (1 – \theta)}{\theta^2} \times \frac{1}{-1}\times (1-\theta)$$
$$ = \frac{1-\theta}{\theta}$$

幾何分布に関しては、べき級数の和を使うことが多いので、毎回書くのも大変なためここで紹介しておく。

$$ 1 + r + r^2 + r^3 + \dots = \sum_{k=0}^{\infty} r^k = \frac{1}{1-r}$$

両辺をrで微分すると、

$$ 1 + 2r + 3\times r^2 + \dots = \sum_{k=1}^{\infty} kr^{k-1} = \frac{1}{(1-r)^2}$$

さらに両辺をrで微分すると、

$$ 2 + 3\times2 r + 4\times 3 r^2 + \dots = \sum_{k=2}^{\infty} (k-1)kr^{k-2} = \frac{2}{(1-r)^3}$$

となる。

幾何分布の分散

$$ E(X(X-1)) = \sum_{x=0}^{\infty} x(x-1) \theta ( 1-\theta)^x $$
$$ = \theta (1-\theta)^2 \sum_{x=2}^{\infty} x(x-1) ( 1-\theta)^x $$
べき級数の和を用いると、
$$ = \theta (1-\theta)^2 \frac{2}{\left[ 1 – (1-\theta) \right]^3 }$$
$$ = \frac{2(1-\theta)^2}{\theta^2} $$

$$V(X)=E(X^2)-[ E(X) ]^2$$
$$ = E[X(X-1)] + E(X) -[ E(X) ]^2$$

$$ = \frac{2(1-\theta)^2}{\theta^2} + \frac{1-\theta}{\theta} – \frac{(1-\theta)^2}{\theta^2} = \frac{(1-\theta)^2 + \theta(1-\theta)}{\theta^2} $$

$$ = \frac{1-\theta}{\theta^2}$$

θの絶対値は1以下なので、幾何分布においては、分散が平均よりも大きくなることがわかる。

7.負の二項分布

r個の不良品が見つかるまでに観測した良品の数Xが従う確率分布で、以下のように表される。

$$P(X=x)={}_{r+x-1} \mathrm{C} _x \theta^r (1-\theta)^x$$
$$ x=0,1,2,\dots $$
負と呼ばれる所以は、以下のように組み合わせを負で表現できるところにある。

以下では実際に負で表現できるかを示す。
$$P(X=x)={}_{-r} \mathrm{C} _x \theta^r(\theta – 1)^x$$

$$ {}_{r+x-1} \mathrm{C} _x $$
$$ = \frac{(r+x-1)(r+x-2)\dots(r+1)r}{x!} $$

分子の数はx個あるので、おのおの-1を掛けると

$$ = (-1)^x \frac{(-r)(-r-1)\dots(-r-(x-2))(-r-(x-1))}{x!}$$

$$ = (-1)^x \frac{(-r)!}{(-r-x)!x!}$$
$$ = (-1)^x {}_{-r} \mathrm{C} _x $$
よって、
$$P(X=x)=(-1)^x {}_{-r} \mathrm{C} _x \theta^r (1-\theta)^x $$
$$ = {}_{-r} \mathrm{C} _x \theta^r (\theta-1)^x $$

負の二項分布の平均

$$E(X) = \sum_{x=0}^{\infty} x \times {}_{r+x-1} \mathrm{C} _x \theta^r (1-\theta)^x $$

$$ =\sum_{x=1}^{\infty} \frac{(r+x-1)!\theta^r (1-\theta)^x}{(x-1)!(r-1)!} $$
ここでx = y + 1 として
$$ =\sum_{y=0}^{\infty} \frac{(r+y)!\theta^r (1-\theta)^{y+1}}{(y+1-1)!(r-1)!} $$
$$ = \frac{r(1-\theta)}{\theta} \sum_{y=0}^{\infty} \frac{(r+y)! \theta^{r+1} (1-\theta)^{y}}{y!r!} $$

$$ = \frac{r(1-\theta)}{\theta} \times $$
$$ \sum_{y=0}^{\infty} \frac{(r+y)(r+y-1)\dots (r+1)}{y!} $$
$$ \times \theta^{r+1} (1-\theta)^{y} $$

$$ = \frac{r(1-\theta)}{\theta} \sum_{y=0}^{\infty} {}_{r+y} \mathrm{C} _{y} \theta^{r+1} (1-\theta)^y $$

ここで負の二項分布の和は1であることから、
$$ = \frac{r(1-\theta)}{\theta} $$
となる。

負の二項分布の分散

$$ E(X(X-1)) = \sum_{x=0}^{\infty} x(x-1) $$
$$ \times {}_{r+x-1} \mathrm{C} _x \theta^r (1-\theta)^x $$

$$ = \sum_{x=2}^{\infty} \frac{(x+r-1)(x+r-2)\dots r}{(x-2)!} \theta^r (1-\theta)^x $$

ここでx = y + 2 として

$$ = \sum_{y=0}^{\infty} \frac{(y+r+1)(y+r)\dots r}{y!} \theta^r (1-\theta)^{y+2} $$
$$ = \frac{r(r+1)(1-\theta)^2}{\theta^2} \sum_{y=0}^{\infty} \frac{(y+r+1)(y+r)\dots (r+2)}{y!} \theta^{r+2} (1-\theta)^{y} $$

$$ = \frac{r(r+1)(1-\theta)^2}{\theta^2} $$
$$ \times \sum_{y=0}^{\infty} {}_{r+y+1} \mathrm{C} _{y} \theta^{r+2} (1-\theta)^y $$

ここで負の二項分布の和は1であることから、
$$ = \frac{r(r+1)(1-\theta)^2}{\theta^2} $$

$$V(X)=E(X^2)-[ E(X) ]^2 $$
$$ = E[X(X-1)] + E(X) -[ E(X) ]^2$$

$$ = \frac{r(r+1)(1-\theta)^2}{\theta^2} + \frac{r(1-\theta)}{\theta} – \frac{r^2(1-\theta)^2}{\theta^2} $$

$$ = \frac{(1-\theta)^2r + (1-\theta)\theta r}{\theta^2} $$
$$ = \frac{(1-\theta)\left[ (1-\theta)r + \theta r \right]}{\theta^2} $$
$$ = \frac{r(1-\theta)}{\theta^2} $$

θの絶対値は1以下なので、負の二項分布においては、分散が平均よりも大きくなることがわかる。

おわりに

以上、これだけ導出を丁寧に書けば後で思い出せるだろうなというレベルで残したつもりですが、冗長的なところもあると思います。おかしなところがあったらご指摘ください。次は連続分布を扱う予定ですが、結構LaTeX書くの疲れますね。

参考文献

[1] 鈴木・山田 (1996), 『数理統計学―基礎から学ぶデータ解析』, 内田老鶴圃
[2] 緑川章一, “負の二項分布”
[3] abrahamcow (2014),”微分方程式によるポアソン分布の導出”
[4] 物理のかぎしっぽ, “定数係数1階線形微分方程式”
[5] 高校数学の美しい物語, “超幾何分布の意味と期待値の計算”

[Stan]ロジスティック回帰の階層ベイズモデルとk-foldsクロスバリデーション

はじめに

stanは意思決定のための分析などでのパラメータ推定に使うことが多く、機械学習のために扱うことはありませんでした。ただ、過去にリク面などでお話したデータサイエンティストの方はstanで機械学習していて、クロスバリデーションもしているとの発言をされていました。
先日、記事を漁っていたらstanでクロスバリデーションを行うためのコードを書いている方を見つけたので、その方のコードをもとに大人のirisであるwineデータを用いて、質の高いワインかどうかを分類するために階層ベイズでのロジスティック回帰モデルを回してみたいと思います。

データについて

UCI Machine Learning Repositoryにある、赤ワインの評価と成分のデータです。データに関する説明はワインの味(美味しさのグレード)は予測できるか?(1)で丁寧になされていますので、ご確認ください。今回は6点以上であれば1を、そうでなければ0を取るものを教師データとしています。

分析方針

今回は階層ベイズモデルを扱うことから、グループごとにロジスティック回帰のパラメータが異なるという仮定を置きます。そのために、citric.acidというデータを3つのカテゴリデータに変換して、それをグループとして扱います。モデルでは、今回のデータセットの変数を全て回帰項として使います。もちろん、回帰用の式からはcitric.acidは除外します。
全体の80%を訓練データに、20%をテストデータとして、10foldsクロスバリデーションでのstanによる予測結果の平均AUCを評価指標とします。最後に、テストデータを用いた予測のAUCを確かめます。また、階層ベイズモデルではないモデルでの10foldsクロスバリデーションでのAUCとも比較します

分析概要

・データ整形
・訓練データとテストデータの分割
・クロスバリデーション用のデータの作成
・stanの実行
・クロスバリデーション結果の出力
・テストデータでの予測
・非階層モデルとの比較

全体のコード以下のリンクにあります。
kick_logistic_regression_allowing_k_hold_cross_validation_hierachical.R
logistic_regression_allowing_k_fold_cross_validation_hierachical.stan

データ整形

階層ベイズで扱うグループをcitric.acidから作っています。

訓練データとテストデータの分割

ワインの質に関するバイナリーデータをこちらで作成し、80%を訓練データに、20%をテストデータに分割しています。

クロスバリデーション用のデータの作成

こちらのコードでは任意の数でクロスバリデーション用のデータを作成し、stanで扱う訓練用データのlistに追加しています。
また、参考にしているブログより転用したstan_kfoldという関数を定義しています。k分割した際のstanの推定結果をリストに格納するための関数です。

stanの実行

こちらのstanのコードでは、M個のグループごとにパラメータが異なるというモデルを書いています。modelブロックの途中でholdoutを入れることで一部のデータを推定に使わないようにしています。

こちらはstanをキックするためのコードです。いつもと違い、先程定義したstan_kfoldを用いています。

クロスバリデーション結果の出力

以下は、k個ずつ手に入ったクロスバリデーションでの推定結果から、各パラメータの平均値を計算し、ロジスティック回帰モデルで2値の予測を行い、平均AUCを計算するコードです。

平均AUCは0.675となりました。すごくいいわけではないですが、手抜きモデルとしてはまずまずと言ったところでしょうか。

テストデータでの予測

以下のコードで最初に分けていたテストデータでの予測結果を返します。

実行の結果、AUCは0.665と、クロスバリデーションでの平均AUCと比べてあまり下がりませんでした。

非階層モデルとの比較

非階層モデルでも同様に10foldsクロスバリデーションの平均AUCを計算しました。非階層モデルよりもAUCが1%ポイントくらいは高いようです。

おわりに

現時点において、stanでの柔軟なモデリングを機械学習に活かす作法について紹介されている文献はあまりなく、選手人口はどれくらいいるのか気になるところですが、発見したブログのやり方でクロスバリデーションをカジュアルに行えるので、より多くの方がstanでの機械学習にチャレンジしうるものだなと思いました。ただ、このレベルの階層ベイズだとrstanarmで簡単にできてしまうので、より深く分析してモデルをカスタムしていきたいですね。

参考情報

[1]Lionel Hertzog (2018), “K-fold cross-validation in Stan,datascienceplus.com”
[2]Alex Pavlakis (2018), “Making Predictions from Stan models in R”, Medium
[3]Richard McElreath (2016), “Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)”, Chapman and Hall/CRC
[4]松浦 健太郎 (2016), 『StanとRでベイズ統計モデリング (Wonderful R)』, 共立出版
[5]馬場真哉 (2019), 『実践Data Scienceシリーズ RとStanではじめる ベイズ統計モデリングによるデータ分析入門』, 講談社

Causal Inference in Economics and Marketingを(今更)読んだ感想と備忘録

2016年にバリアンがPNAS(Proceedings of the National Academy of Sciences of the United States of America)に投稿した”Causal Inference in Economics and Marketing”というペーパーを見つけました。内容的に非常にわかりやすかったし学びがあったので、要約して自分のためにメモを残しておきたいと思います。要約なんぞいらないという方はこちらのPDFを見ていただけると良いです。

修士が終わってからWebマーケティングの仕事をずっとしていますが、ABテストができるケースが多かったので因果推論が必要になる場面はあまりありませんでした。ただ、事業会社でしばしば行われるオフライン広告や社内の研修の効果に関してはABテストのような実験はできませんので、因果推論のニーズは「実は」高いと思います。ただ、広告を打つ主体・研修を進める主体がその分析のニーズがあるかというと全然そんなことはないと思います。今回のバリアンの素晴らしい記事をもとに因果推論の民主化をできると、説明の手間が省けるぶん世のデータサイエンティストは少しだけ楽になるかもしれません。

目次

・欠落変数の問題
・因果推論の重要なコンセプト
・因果効果を推定できる方法
 ・自然実験
 ・操作変数法
 ・回帰分断デザイン
 ・Difference in differences(DID)
・おわりに
・参考文献

欠落変数の問題

まずは以下のようなシンプルな線形モデルについて振り返ってみます。y_cはある都市でのサーフィン映画の一人あたりの売上、x_cはある都市でのサーフィン映画に関する広告支出を、bは係数、e_cは誤差項を表します。

この式を推定できれば広告主は何が嬉しいのかと言うと、1単位あたり広告支出を増やせば、どれくらいサーフィン映画の売上が増えるのかがわかるということです。
しかしながら、様々な都市を混ぜこぜしたデータからこの式を推定したものを使って、求まったbが10だとします。その際に、100円突っ込めば1000円儲かる、濡れ手で粟とはこのことや!状態が訪れるのでしょうか。
回帰分析には欠落変数の問題があり、欠落変数が説明変数と相関していると、濡れ手で粟とは言えない状況が起きうります。

まず、本来モデルに加えるべき変数がモデルに含まれていない場合、その変数はすべて誤差項に含まれることになります。
係数bの数式を確認してみると、

から、yを代入すると、

となります。右辺のbは真のbの係数とします。もし、誤差項eと説明変数が相関していれば、右辺二項目の影響により偏った推定量を得ることになります。つまり、本来モデルに加えるべき変数が漏れており、それが説明変数と相関している場合は、濡れ手で粟状態になるとは言えず、100円広告に突っ込んでも1000円儲かると期待することはできなくなります。せっかく推定したのに大損するという、非常に気の毒な結末があるでしょう。

釈迦に説法ですが、被説明変数と説明変数の両方に対して相関してしまう変数を交絡変数と呼びます。今回のサーフィン映画の場合、モデルに含まれていない変数として「サーフィンへの関心度」というものが考えられます。例えば、ハワイ州のホノルルのようなエリアではサーフィンが盛んですしサーフィン映画も人気になりえます。他方、ノースダコタ州のファーゴのような、群馬・栃木・埼玉・岐阜のような海なし県のようなエリアの場合、サーフィンもサーフィン映画も人気ではないでしょう。

ここでのサーフィンへの関心度は「サーフィン映画に対する広告支出」にも影響を与える可能性もありますし、そもそもの「サーフィン映画の興行収入」にも影響を与える可能性があります。

また、広告費自体は、映画の広告担当者が様々なドメイン知識を利用して意思決定しているのが普通で、都市に対してランダムに広告予算を割り振っているという、無能な運用者を想定することは難しいです。そのため、都市間でランダムに広告費が割り振られると考えるような分析者はあまり賢明ではないでしょう。現実世界において、人間が意思決定した結果のデータのほとんどが交絡変数を持つと考えられます。

因果推論の重要なコンセプト

因果推論における重要なコンセプトは、反実仮想と実際の結果との比較をすることです。この場合の反実仮想は処置されなかった場合の処置群の結果となります。

ランダム化実験が可能な場合、セレクションバイアスはゼロになります。そのため、ABテストなどでランダム出し分けができるWebマーケティングは因果推論が行いやすいのです。しかしながら、実務においてのランダム化実験や、ランダムでなくともそもそもの実験自体がハードルが高いケースが多いです。以降ではそのような現状を踏まえた上で、因果効果を推定するための方法を紹介します。

因果効果を推定するための方法

自然実験

自然実験は、手元のデータからランダムに近そうな事象があるのではないかを見出すというアプローチです。バリアンの例では、スーパーボウル(アメリカンフットボール)の試合での映画広告出稿に関して、ホームシティのチームの試合がそのスタジアムで開催される場合に、10~15%ほど観客が多いという前提があるとして、映画広告の出稿の意思決定をする際にどのチームがそのスタジアムで試合が行われるかわからないという点に着目していました。つまり、どのチームがどのスタジアムで試合するかがランダムに決まると考え、広告を打たなかった場合の売上という反実仮想を推定することで、セレクションバイアスなく広告による映画の売上げ増の効果を推定することができます。

操作変数法(二段階最小二乗法)


一般的に広告費をどれだけ支払うかの意思決定をしている人は売上に影響を与えるような変数を考慮して、広告費を決めていると考えられます。しかしながら、それが分析者のモデルの中で考慮されず、誤差項として扱われてしまうことがあります。そのような状況下で適切な広告効果を推定する方法として操作変数法があります。再び釈迦に説法ですが、操作変数は広告費を通してのみ売上に影響を与える変数のことを指します。ここで、スーパーボウルの話に戻すと、プレイオフなどでのホームチームの勝利はそのスタジアムでの広告費の増加につながると考えられます。しかしながら、プレイオフに行けるかどうかは映画の売上に影響を与える可能性は低いと考えられます。つまり、プレイオフに進んだかどうかは操作変数として扱うことができそうです。

操作変数法によるbの推定量は以下のようになります。

操作変数zと誤差項eが直交することで最後の項は消え、真のパラメータ、つまり真の広告効果を求めることができます。あるいは、以下のような2つの方程式を用いることもできます。

第1段階ではプレーオフで勝ったかどうかを表すz_cを用いて、2つ目の式のパラメータを推定します。第2段階では第1段階で推定したaを使って予測したx_cを説明変数としてbを推定します。このようなアプローチを字面の通り、二段階最小二乗法と呼びます。この推定結果は操作変数法と等しくなります。

回帰分断デザイン(Regression discontinuity)

閾値を超えることによる因果効果に関心があるケースです。閾値よりも低いもの、高いものの両サイドに近いものを比較する。例えば、学校のクラスで40名のクラスと41名のクラスがいたとして、教育パフォーマンスのクラスの大きさによる因果効果を検証したい場合や、住宅周辺のブロードバンドのスピードが閾値を超えたと閾値を下回るエリアから、ブロードバンドの速度が住宅価格に与える因果効果を検証したい場合などに用いられるます。学校のクラスで40名クラスや41名クラスに配属されるというのはランダムに近いと考えます。


(ペーパーの図を引用しております。この図では21歳前後かどうかという歳のわずかな違いにおいてランダムに近い、処置群とコントロール群が生成されていると考えます。)

Webマーケティングの現場での活用においては、特定のページの閲覧数がある閾値を超えた場合に、そのユーザーに対して特別に処置(例えばクーポンを配る)を行うケースにおいて、閾値を超えたか超えなかったかの比較から処置の因果効果を推定することになります。つまりWebページの閲覧数が閾値を超える超えないはランダムに近いとしています。その際、閾値を超えたものの処置をしそこねた(クーポンを配り損ねた)というのがWeb広告における反実仮想となります。

バリアンのペーパーによると、Google、Microsoft、Facebookのエンジニアはこれらの実装を行い、閾値による比較を行っているようです。

Difference in differences(DID)

Difference in differences(DID)は時系列データでの因果推論において役に立つ手法です。
ここでは介入が仮になかったとしたら、処置群とコントロール群の前後の差は等しいという、以下のような仮定を置きます。

y_TFは処置群が、介入を受けなかった場合の後の結果(つまり反実仮想)、y_TBは処置群が介入を受ける前の結果、y_CAはコントロール群の後の結果、y_CBはコントロール群の前の結果を表しています。
これを処置群の後の結果y_TAから差し引くと以下のようになります。

つまり、処置群の前後の差は、処置群とコントロール群の前後の差が等しいのであれば、処置群の前後の差からコントロール群の前後の差を差し引いたもので表現することができます。これがDifference in differencesです。変動を見たい場合はブートストラップ法を適用することで変動幅を求めることができます。

水準のみならず、前後の比に関しても同様に扱うことができ、その場合は以下のように対数の差となります。

また、DIDにおいては反実仮想を推定する方法として回帰モデルが用いられることがあります。例えば、線形回帰やランダムフォレスト回帰などです。アプローチとしては機械学習に近く、共変量(その日の天気、ニュースなどのイベント、その他の外生変数)を用いてコントロール群の一部のデータで訓練します。続いて、残りのコントロール群のデータを用いてテストをし、モデルとしての汎化性能を高めます。そのモデルを用いて、介入期間における反実仮想を予測し、その反実仮想と実績を比較することで因果効果をもとめます。


(ペーパーの図を引用しております。こちらでは訓練データを用いて推定したモデルでテスト期間を予測し、介入期間での反実仮想を予測により求め、介入後の実績値と比較しています。)

おわりに

因果推論に関して既に様々な文献がありますが、バリアンの説明が一番頭に入ってきやすかった気がします。もちろん、因果推論のすべてのアプローチが網羅されているわけではないですが、これから因果推論の勉強を頑張ろうというモチベーションを上げるための教材としては良いなと思いました。

参考文献

[1] Varian(2016), “Causal Inference in Economics and Marketing”, PNAS
[2] “回帰不連続デザインRegression discontinuity design(および分割時系列デザイン)”
[3] 岩波データサイエンス刊行委員会(2016), 『岩波データサイエンス Vol.3』