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

はじめに

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

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

今回登場する連続分布

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

  • 1.一様分布
  • 2.正規分布
  • 3.対数正規分布
  • 4.ガンマ分布
  • 5.指数分布
  • 6.ワイブル分布
  • 7.ベータ分布
  • 8.コーシー分布

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

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

参考文献

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

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

はじめに

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

※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』

[Python]機械学習などでテキストデータを特徴量にする際のソースコード集

テキストデータの特徴量化について

仕事ではテキストデータを多用するので、機械学習などで扱うためにテキストデータを特徴量にするためのアプローチを色々と整理してソースコードを残しておきたいと思います。今回はあくまでも私の知っているものだけなので、網羅性はないかもしれませんが悪しからず。
(2019/08/18 追記)Stackingをカジュアルに行えるvecstackというモジュールを用いた予測も試してみました。下の方の追記をご覧ください。

アプローチ

テキストデータを特徴量にする際のアプローチとしては、以下の3つが良く使っているものとなります。
・単語ベース
・クラスタ、トピック、分散表現ベース
・文書間の類似度ベース

今回扱うデータ

ひょんなことから、昨年10月くらいに取りためたマンションの施設情報のテキストです。

緑色が印象的な某不動産紹介サイトをクローリングしました。全部で1864件ほどの文書数となります。

加えて、デザイナーズマンションかどうかのフラグを作成しました(17%くらいがデザイナーズマンションの割合)。これでもって、マンションの施設情報からデザイナーズマンションかどうかを分類できるかチャレンジしたいと思います。
ここにデータを置いていますので、興味のある方はご利用ください。

今回扱うモデル

ランダムフォレストです。10foldsクロスバリデーションによるAUCの結果を各手法のスコアとして扱います。

こちらは、任意の手法に関して10foldsクロスバリデーションで実行し、AUCのグラフを生成してくれるソースコードです。主にscikit-learnのサイトに載っているものです。引数のclassifierをsklearnの任意のモデルのインスタンスで渡せば動きます。

単語ベース

シンプルに単語をそのまま特徴量にするというものですが、文書によっては単語数が多すぎて収集がつかないと思います。そこで単語を簡単に選択できるDocumentFeatureSelectionというパッケージを利用します。

このパッケージでは
・TF-IDFベースの特徴量選択
・PMI(Pointwise Mutual Information)ベースの特徴量選択
・SOA(Strength of association)ベースの特徴量選択
・BNS(Bi-Normal Separation)ベースの特徴量選択
を行うことができます。

まずは今回のベースラインとして、単語のカウントベースでの特徴量を扱いたいと思います。
その前に、GitHubに上がっているデータに対して以下のように簡単な前処理をしておきます。

ようやくベースラインの予測となります。以下のコードを実行すると、ROCが描かれた図がJupyter上で表示されます。

AUC82%というのはベースラインとしてはなかなか強敵なのではないでしょうか。

さて、本題の特徴量選択パッケージの適用をするためのソースコードを以下に記します。

以上のソースコードを実行すれば、tf_idf_scored_df、pmi_scored_df、soa_scored_df、bns_scored_dfにスコアを付与された単語のリストが手に入ります。

ここでは各スコアに関してアドホックに閾値を設けて、特徴量として利用することにします。

TF-IDFベースの特徴量選択

PMIベースの特徴量選択

SOAベースの特徴量選択

BNSベースの特徴量選択

クラスタ、トピック、分散表現ベース

続いて、k-meansやLDAやword2vecを用いて特徴量を作成する方法です。今回はk-means、ミニバッチk-means、LDA、FastTextによる分散表現を扱います。

k-means、ミニバッチk-means

LDA

こちらはgensimでLDAを推定し、推定したトピックの割合をデータフレームで返すコードです。

トピック数をとりあえず30個くらいに指定して推定したトピックの割合を特徴量として文書分類を行います。そのため、特徴量の数は30個になります。

FastTextによる分散表現

今回はデータ数が少ないことから、学習済みの分散表現を用います。日本語のコーパスに対して、FastTextで推定された分散表現となります。学習済み分散表現はこちらから拝借しました。

分散表現は単語に対して計算されるので、単語に対して分散表現を足し合わせたものを特徴量として扱います。ここでは分散表現の合計値、平均値、TF-IDFで重みを付けた平均値の3つのパターンを試します。

合計値ベース

平均値ベース

TF-IDFで単語を重みづけた平均値ベース

文書間の類似度ベース

今回は、デザイナーズマンションの定義文に似ているかどうかという観点で類似度ベースの特徴量を作ってみたいと思います。

今回は変数が一つだけなので、機械学習はせず、デザイナーズマンション割合との関係を図示するにとどめておきます。横軸がデザイナーズマンションの定義と施設情報の類似度で、縦軸がデザイナーズマンション割合です。

どうやら、途中でデザイナーズマンション割合がピークを迎えるようです。

おわりに

最先端の手法は調べれていないですが、テキストデータを特徴量に落とし込む手段を備忘録として残しておきました。今回あげた中では、SOAベースの特徴量選択のAUCが83%と一番高かったですが、ベースラインが82%と僅差でした。そして、分散表現形のものは80%に届いた程度です。余力があれば新しい特徴量の作り方が分かり次第アップデートしようと思います。

追記

“Automate Stacking In Python How to Boost Your Performance While Saving Time”という記事を見つけたので、紹介されているvecstackモジュールを使って今回のモデルに関して簡単にstackingしてみようと思います。
コードに関しては、こちらのGitHubに上げています。試してみた所、AUCは88%になりました。結構上がりましたね。しかもコードはめちゃ短いので楽です。

参考文献

[1]Julian Avila et al(2019), 『Python機械学習ライブラリ scikit-learn活用レシピ80+』, impress top gear
[2]Receiver Operating Characteristic (ROC) with cross validation
[3]@Kensuke-Mitsuzawa(2016), “テキストデータで特徴量選択をする”, Qiita
[4]JapaneseTokenizer 1.6
[5]DocumentFeatureSelection 1.5
[6]自然言語処理における自己相互情報量 (Pointwise Mutual Information, PMI)
[7]【Techの道も一歩から】第3回「第11回テキストアナリティクス・シンポジウム」
[8]文書分類タスクでよく利用されるfeature selection

ベイジアン線形回帰モデルの式変形とRでのギブスサンプリングの適用

今回は特に目新しい手法というわけでもなく、線形回帰モデルのギブスサンプリングについて忘備録として残しておきたいと思います。
ベイジアン線形回帰モデルはプログラミング言語で言う、Hello World!的なものなので、あえてブログで取り上げる必要があると考えていないのですが、導出をしては忘れの繰り返しが嫌なので自分のために残しておこうと考えました。加えて、Stanのありがたみを感じられ、Stanへのコミットメントが増すのではないかとも期待しています。

・モデル
・数式の展開
・Rのコードの紹介
・おわりに
・参考情報

モデル

東北大学の照井教授の『ベイズモデリングによるマーケティング分析』に載せられている表記に従い、以下のように記します。

説明変数の数がk個の正規線形モデル

を考える。その場合、尤度関数は

となる。

係数パラメータの事前分布や条件付きの誤差分散の事前分布は以下のように設定する。(βは正規分布に従い、σ2|βは逆ガンマ分布に従う。)

数式の展開

私が大学院生だった時に、数式の展開をどう進めるかを手っ取り早く知る方法としては、「ネットに上がっている海外の大学院の講義資料を漁る」という作戦を取っていました。こうすることで数学のセンスがそれほど高くなくても、理解し進めることができました。今回に関してもおそらく、わかりやすく解説している海外の研究者がいるはずだと思い漁ってみたところ、コロンビア大学の機械学習の講義資料を見つけることができました。

資料はこちらのPDF(Course Notes for Bayesian Models for Machine Learning)で126ページにもなっていますが、導出のステップなどが非常に丁寧に書かれています。

それでは、今回の講義ノートを参考にしながら、線形モデルにおいて、ギブスサンプリングを行うところまでの導出を行いたいと思います。

まず、同時事後分布を以下の左辺のように置き、ベイズの定理を用いて右辺のように表記する。

次に、条件付き確率の定義と先程の尤度関数から以下のようになる。

yが与えられたもとでのp(y)は一定のため、比例している分子だけを残すと以下のようになる。

同時事後分布に事前分布の関数を代入していくと、

となる。両辺について対数を取ると、

となる。ここでβやσ2についての事前分布の形状から、同時事後分布におけるβやσ2について整理するための目標となる形状を確かめる。
まず、βはp(β)の定義より、対数を取りβについて整理すると、

となる。つまり、1/B0や1/B0・β0に該当する表現を先程の対数を取った同時事後分布から得ることを目標とする。
他方、σ2についても同様に、p(β|σ2)の定義より対数を取りσ2について整理すると、

となる。つまり、ν0やδ0に該当する表現を、同じく対数を取った同時事後分布から得ることを目標とする。

以上のパラメータごとの目標とする形状になるように各々のパラメータについて、対数を取った同時事後分布を整理する。

まずはβについてまとめ、関係のない項をconst.にする。

先程もとめた目標の形状を当てはめると以下のようになる。

よって、βの事後分布は以下のようになる。

他方、σ2についても同様に、関係のない項をconst.にし、目標の形状にまとめると以下のようになる。

目標の形状と比較すると以下のようになる。

よって、σ2の事後分布は以下のようになる。

Rのコードの紹介

条件付き事後分布からβやσ2の従う分布の形状がわかったので、それらを使ってRでギブスサンプリングを行います。先日、たまたま見つけた線形回帰モデルのギブスサンプリングのRのソースコードを拝借しようと思います。

ギブスサンプリングでは、先程導出した条件付き分布からβ→σ2と交互にサンプリングしていきます。それを記述したRコードは以下の通りです。

先程導出したβの事後分布である正規分布からのサンプリングの後(15~17行目)、そのサンプリングしたβを用いて、同じく先程導出したσ2の事後分布である逆ガンマ分布からサンプリングし(19~21行目)、それを指定した回数だけ繰り返し、所定の数まではバーンインとして除外します。(25~27行目)こうして導出した数式と、ギブスサンプリングのコードを見比べると理解が捗ると思いました。

実際に、先程のGitHubのソースコードを回してみると、以下のようにギブスサンプリングのイタレーションのプロットや、パラメータの事後分布を確認できます。

全体のコードはこちらです。

おわりに

シンプルなモデルですらこれだけ導出に手間がかかるということからも、Stanなどの確率的プログラミング言語のありがたみは非常に大きいなと思いました。こうして残すことで今後忘れたとしてもすぐに思い出せる気がします。
しかしながら、Stanでは事前分布と尤度を指定してしまえば、事後分布を計算し、知りたいパラメータについて解いた条件付き分布からサンプリングしてくれるわけですから、研究者の寿命を伸ばしたと言っても過言ではないかもしれません。

参考情報

[1]John Paisley (2016), “Course Notes for Bayesian Models for Machine Learning”, Columbia University
[2]照井伸彦 (2008), 『ベイズモデリングによるマーケティング分析』, 東京電機大学出版局
[3]須山敦志 (2017), 『機械学習スタートアップシリーズ ベイズ推論による機械学習入門』, 講談社
[4]stablemarkets,BayesianTutorials/MultipleLinearReg/multiplelinearreg.r

[Stan]生存時間分析のコードと便利なデータセットについて

はじめに

仕事で生存時間分析を使うことは結構あるのですが、マーケティングの良いデータセットがない印象でブログにしにくいと感じていました。また、Stanでの生存時間分析の事例もあまり把握していません。そこで使えそうなデータセットやStanのコードを探して、そのデータに対して生存時間分析を適用してみたいと思います。

目次
・生存時間分析とは
・生存時間分析で使えるデータ
・生存時間分析をマーケティングで使う際の用途
・先行研究
・生存時間分析で使えるデータセット
・Stanでの実行例
・おわりに
・参考文献

生存時間分析とは

生存時間分析は、ある時点から興味のあるイベント(マーケティングだと解約など)が発生するまでの期間を分析対象としています。データを手に入れた時点で、すでに解約して、真の累積の契約期間が判明している人と、解約しておらず今後いつ解約するかわからない中での累積の契約期間が残されている人のようなデータを扱うことが多いです。ここでの後者をcensoring(打ち切り)されたデータと呼びます。

生存時間分析をマーケティングで使う際の用途

ブログなどを読み漁る限り、以下の用途で生存時間分析を活用できるようです。

  • 顧客のサービス離脱率の予測、離脱原因の特定
  • 顧客がマーケティングキャンペーンに反応するまでの期間の長さ
  • 故障率の予測、故障原因の特定

先行研究

Stanを用いた分析事例は、調べた限りですが以下のモデルがあるようです。

  • 指数分布のモデル
  • Weibull(ワイブル)分布による比例ハザードモデル
  • ハザードの対数値についてのランダムウォークモデル
  • 2階差分のマルコフ場モデル(生存時間の確率分布は正規分布)
  • 1階差分のランダムウォークモデル(生存時間の確率分布は正規分布)

生存時間分析で使えるデータセット

事例を調べる過程で見つけた、生存時間分析に適したデータセットは以下の通りです。

どうやら、マーケティング、HR、離婚、再犯と幅広いデータがオープンソースで手に入るようです。

Stanでの実行例

今回は、「Using Customer Behavior Data to Improve Customer Retention」のデータセットを用いて、先行研究のソースコードにより生存時間分析をしてみようと思います。データは電話会社の顧客の解約に関するもので、様々な顧客の履歴データなどが用意されています。
先行研究のソースコードはWeibull分布を想定した比例ハザードモデルです。今回は決済の電子化の有無と離脱の関係を確かめてみます。なお、今回の打ち切りデータは契約期間となります。

まずはStanのコードはこちらです。Xobs_bgに説明変数が来るようにデータを用意しておきます。

続いて、このStanコードをキックするためのRのソースコードです。 元のデータが7043件と多すぎるのでランダムサンプリングしています。サンプリング数を8000、チェイン数を4にして実行します。(なお、可視化のソースコードもあるので結構長くなっていますので。頑張ってスクロールしてください。)

Rhatは全て1.05以下になっています。

traceplotを見る限り、重なり合っているので問題なさそうです。

各パラメータごとの自己相関係数です。こちらも問題なさそうです。

推定したパラメータの分布です。

横軸は推定した継続期間です。決済の電子化をしていない消費者は、契約期間の短い際の確率密度が低い傾向があるようです。

どうやら離脱率に関して決済の電子化をしていない消費者は、そうでない消費者よりも低いようです。

こちらは95%で取りうる範囲をそれぞれプロットしたものです。

おわりに

Stanで生存時間分析を行うという事例はそんなに多くはないものの、業界の長たちが良いコードを作成してくれていました。また、面白そうなデータセットも見つけることができました。このようなデータがもっと広まっていけば、マーケティングにおける生存時間分析がより活発に行われるのかもしれません。

参考文献

[1] 豊田秀樹 (2017) 『実践 ベイズモデリング -解析技法と認知モデル-』 朝倉書店
[2]生存時間解析入門
[3]比例ハザードモデルはとってもtricky!
[4]Stanで生存時間解析(Weibull 回帰)
[5]生存時間分析をStanで実行してみた
[6]階層ベイズ生存解析を用いたwebサイトの訪問者分析に関するStanでの実装
[7]生存時間分析 – ハザード関数に時間相関の制約を入れる
[8]Bayesian Survival Analysis 1: Weibull Model with Stan
[9]Bayesian Inference With Stan ~062~
[10]生存時間解析について – 概要編
[11]Survival Analysis for Employee Attrition ※kaggleで提供されているHR系のデータをサバイバル分析に用いている。
[12]Survival Analysis with R※Random Forests Modelによる生存時間の推定がなされている。
[13]Survival Analysis with R and Aster ※服役後の犯罪に関する分析や、離婚の分析などをしている。
[14]Survival Analysis of Mobile Prepaid Customers Using the Weibull Distribution(ダウンロード注意)

[Stan]項目反応理論(IRT)の段階反応モデルでbaysemのアンケートデータの分析をしてみる

はじめに

stanのユーザーガイドを見ていて、項目反応理論(IRT)についての章があり気になりました。勉強会のLTなどで手法の名前をちらっと聞いたことはあったのですが、使い道について調べていませんでした。ビジネスにおける実活用もしやすそうだと思ったので、カジュアルに分析して備忘録として残したいと思います。

目次
・項目反応理論(Item Response Theory:IRT)とは
・ビジネスでの適用可能性について
・データ
・モデルの推定
・結果の解釈
・おわりに

項目反応理論(Item Response Theory:IRT)とは

関西学院大学の教授のブログによると、

項目反応理論とは、テストについての計量モデルで、問題に対する正解・不正解のデータから、問題の特性や、回答者の学力を推定するためのモデルです。

とあります。また、Wikipediaによると、TOEFLの問題の評価のために使われているそうです。

主に、バイナリーと順序変数のモデルがあるようで、以下の母数がモデルに想定されています。どちらもほぼ同じです。

回答が2値変数のモデル

  • 2母数のロジスティックモデル
    • 特性値(例えば、広告配信の満足度とか)
    • 識別度母数(項目特性曲線の傾き)
    • 困難度母数(項目特性曲線の切片)
    • 定数

回答が順序変数のモデル(まずい < まぁまぁ < おいしい)

  • 段階反応モデル
    • 特性値(例えば、広告配信の満足度とか)
    • 識別度母数(項目特性曲線の傾き)
    • 困難度母数(項目特性曲線の切片)

※項目特性曲線は横軸に特性値、縦軸に質問の正答率を取ったものです。

ビジネスでの適用可能性について

  • 顧客のアンケート結果の解釈
    • 異質な集団間の得点を比較可能
    • 異なる尺度間の得点を比較可能(昔のアンケートだと5段階、今のアンケートは7段階などの状況はビジネスデータでありうる。)
  • 人事評価のバイアスの統制
    • 採用面接時の個人特性の正当な評価
  • アンケート項目の項目削減によるアンケートコストの低減
    • 各アンケート項目が理解されたかどうかを分析し、一つ一つのアンケート項目の精度を高める

データ

今回扱うデータはbaysemパッケージに入っているデータセットです。Yellow Pagesの広告プロダクトにおける満足度サーベイの回答データで、全ての回答は1から10のスケールで点数が付けられています(1がPoorで10がExcellent)。質問数は10個で、回答数は1811件です。

各質問の内容(baysemパッケージのドキュメントに載っていました。)
q1:全体の満足度

価格について
q2:競争的な価格設定
q3:昨年と同じ広告の最小値に対しての価格の引き上げ
q4:消費者の数に対しての適切な価格設定

効果について
q5:広告の購入の潜在的な影響
q6:広告を通じて自身のビジネスへの集客ができたか
q7:多くの消費者にリーチしたかどうか
q8:年間を通じて消費者に対する長期での露出があったか
q9:多くの家計やビジネスを必要としている人に届いたかどうか
q10:ビジネスを必要としている地理上のエリアに届いたかどうか

今回のIRT適用における特性値は、「広告プロダクトに関する満足度の傾向」としてみたいと思います。

モデルの推定

今回は教科書にならって以下の段階反応モデルを用います。

ここでaは識別力(広告の満足度が高まりやすいかどうか)、bは境界パラメータ(回答カテゴリ間の境界値)、θは特性(回答者がどれだけ広告に満足しているか)を表しています。Dは定数項で、以下では1とおいています。cはアンケートの回答のカテゴリ番号です。今回の例では10段階の評価が入ることになります。最後に、uは反応を、jは質問の番号を表しています。

実践 ベイズモデリング -解析技法と認知モデル-

こちらの本のサポートサイトからダウンロードできるzipファイルにstanのコードやRコードがありますので、そちらを利用しています。

モデルですが、以下のような設定となっています。

こちらをキックするためのRコードです。

処理時間としては、2014年末モデルのMacbook Proのcorei5、メモリ8GBで数時間程度でした。(正確な時間はわかりませんが、寝て起きたら計算が終わっていました。)


どうやら収束してそうです。


Rhatも1.1未満におさまっています。

結果の解釈

まず、推定した特性値の値のユーザーごとの平均値を求めて、ヒストグラムを描いてみます。どうやら、上限周辺にやたらと高い評価をしてそうなユーザーがいるようです。

最後に、項目特性曲線を質問ごとに、そして回答ごとに描いてみようと思います。

質問1~10に関して、10段階の回答ごとの項目反応曲線を以下に描いています。上まで戻るのが面倒なので、質問内容を再掲します。

q1:全体の満足度
q2:競争的な価格設定
q3:昨年と同じ広告の最小値に対しての価格の引き上げ
q4:消費者の数に対しての適切な価格設定
q5:広告の購入の潜在的な影響
q6:広告を通じて自身のビジネスへの集客ができたか
q7:多くの消費者にリーチしたかどうか
q8:年間を通じて消費者に対する長期での露出があったか
q9:多くの家計やビジネスを必要としている人に届いたかどうか
q10:ビジネスを必要としている地理上のエリアに届いたかどうか

これらの傾向から、9〜10点を獲得するにはある程度は特性値が高まる必要がある質問としては、q1〜q6のように見えます。価格や購買など自身のビジネスに直結しそうな質問が多い印象です。逆にふわっとした質問であるq7~q10は特性値が低くても9〜10点を取れる可能性が高い傾向があります。

おわりに

Stanのユーザーガイドを読むことで、普段自分が業務で扱っているアプローチなどが如何に限定的であることが実感できました。今回はIRTのアンケートデータへの適用事例を知れ、そこから様々な文献や便利なコードに至ることができました。社内のアンケートデータへの適用は面白そうだと思いますので業務で使ってみようと思います。

参考情報

[1] 豊田秀樹 (2017) 『実践 ベイズモデリング -解析技法と認知モデル-』 朝倉書店
[2] Yoshitake Takebayashi (2015) 「項目反応理論による尺度運用」 SlideShare
[3] 持主弓子・今城志保 (2011) 「IRTの組織サーベイへの応用」
[4] 清水裕士 (2017) 「項目反応理論をStanで実行するときのあれこれ」 Sunny side up!
[5] 清水裕士 (2016) 「Stanで多次元項目反応理論」 Sunny side up!
[6] 小杉考司 (2013) 「項目反応理論について」
[7] Daniel C. Furr et al. (2016) “Two-Parameter Logistic Item Response Model”
[8] daiki hojo (2018) “Bayesian Sushistical Modeling” Tokyo.R#70
[9] abrahamcow (2017) 「[RStan]項目反応理論の応用でフリースタイルダンジョン登場ラッパーの強さをランキングしてみた」

データアナリストがBIダッシュボードのお手伝いをする前の調べ物

はじめに

これまで、データ分析・機械学習などの知識を重点的に学んできたので、BI周りはかなり疎かになっていました。今の知識のままではBI案件のお手伝いが大変だなと思ったので、少しはリソースを割いて勉強していこうと思います。今回はコードも何も出てこない内容なので、データサイエンティストの方はそっ閉じで良いと思います。
海外のブログで、“12 Best Business Intelligence Books To Get You Off the Ground With BI”という記事があり、”Performance Dashboards: Measuring, Monitoring, and Managing Your Business (English Edition)“という書籍が紹介されていました。それについてわかりやすくまとめられたスライドがあったので、2006年と古いですが、何もBIを知らないものとしては非常に勉強になると思い、訳して自分のための忘備録としておきます。

目次

・パフォーマンスダッシュボード
・ダッシュボードの不満
・パフォーマンスダッシュボードの2つの原則
・戦術的なドライバー
・戦略的なドライバー
・パフォーマンスダッシュボードは何で構成されるか
・おわりに

パフォーマンスダッシュボード

  • ダッシュボード
  • パフォーマンスチャート
    からなる

ダッシュボードの不満

意思決定を阻害するもの

  • データが多すぎる
  • 情報が少なすぎる
  • 届くのが遅すぎる

パフォーマンスダッシュボードの2つの原則

パフォーマンスダッシュボード = BI + 企業の経営マネジメント

(画像はPDFより拝借しております。)

  • 1.BI
    • 情報(DWH)
    • 知識(分析ツール)
    • 計画(ルール、モデル)
    • 行動(レビュー、計測、洗練)
    • 知恵
    • イベント発生→情報へ
      の繰り返し
  • 2.企業の経営マネジメント
    • 戦略(ミッション、バリュー、ゴール、目的、インセンティブ、戦略マップ)
    • 計画(予算、計画、見込み、モデル、イニシアティブ、ターゲット)
    • モニタリング、分析(パフォーマンスダッシュボード)
    • 行動、調整(行動、決定、見直し)

戦術的なドライバー

  • 利用者と共鳴する
    • 一つのスクリーンでいくつかの領域のステータスをモニタリングできる
    • 重要な指標のグラフ表示
    • 例外的な状況に関してアラートを上げる
    • クリックして分析し、詳細を深掘りできる
    • ルールに基づきカスタマイズされた表示
    • 訓練が要求されない
  • リッチなデータ
    • 複数の情報ソースからブレンドされたデータ
    • 詳細も集計値もある
    • 履歴もリアルタイムのデータもある
  • 労働者に力を与える
    • 本当に重要なことにユーザーを集中させる
    • 労働者の貢献がどのように集計されているかを示す
    • ゴール、競争、インセンティブで動機づけをする
    • プロアクティブな介入を促進する

戦略的なドライバー

  • ビジネスを調整する
    • 皆同じデータを使う
    • 皆同じ指標を使う
    • みな同じ戦略で働く
  • コミュニケーションの改善
    • コミュニケーション戦略のためのツール
    • マネジャーとスタッフのコラボレーション
    • 部門間のコーディネート
  • 視認性とコンプライアンスの向上
    • 驚きの少なさ
  • 戦略的なドライバーの5つのC
    • Communicate
    • Compare
    • Collaborate
    • Coordinate
    • Congratulate

パフォーマンスダッシュボードは何で構成されるか

  • 3つのアプリケーション
  • 情報の3つのレイヤー
  • パフォーマンスダッシュボードの3つのタイプ

3つのアプリケーション

  • モニタリング
  • 分析 
  • コラボレーション


(画像はPDFより拝借しております。)

情報の3つのレイヤー


(画像はPDFより拝借しております。)

  • モニタリング:グラフ、図形、チャート
  • 分析:ディメンション、階層、細かく分割
  • レポーティング:DWHのクエリ実行、運用レポート
  • これらをプランニングする
     ・計画、モデル、予測、更新

ダッシュボード

  ・目的:現在の活動状況を測る
  ・ユーザー:経営者層、マネジャー、スタッフ
  ・更新頻度:即時
  ・データ:イベント
  ・クエリ:リモートシステムで実行
  ・画面:チャート

スコアカード

  ・目的:進行状況を示す
  ・ユーザー:経営者層、マネジャー、スタッフ
  ・更新頻度:周期的なスナップショット
  ・データ:サマリー
  ・クエリ:ローカル環境のデータマートで実行
  ・画面:図形

・経験則
 ・ビジネスユーザーが好むものは何でも使う!

パフォーマンスダッシュボードの3つのタイプ

  • 業務系
    • 焦点:モニタリング業務
    • 重点:モニタリング
    • ユーザー:管理者
    • スコープ:現場
    • 情報:詳細
    • 更新頻度:日中
    • 適しているのは:ダッシュボード
  • 戦術系
    • 焦点:プロセスの最適化
    • 重点:分析
    • ユーザー:マネジャー
    • スコープ:部門
    • 情報:詳細/サマリー
    • 更新頻度:日次/週次
    • 適しているのは:BIポータル
  • 戦略系
    • 焦点:戦略実行
    • 重点:コラボレーション
    • ユーザー:経営者層
    • スコープ:企業
    • 情報:サマリー
    • 更新頻度:月次/四半期
    • 適しているのは:スコアカード

パフォーマンスダッシュボードをどのように作るか?

3つのアーキテクチャー

・ビジネスアーキテクチャーとテクニカルアーキテクチャー
・BIアーキテクチャー
・データアーキテクチャー

ビジネスアーキテクチャー

  • ステークホルダー:投資家、取締役、全従業員、顧客、サプライヤー、監督機関
  • 戦略:ミッション、ビジョン、バリュー、ゴール、目的、戦略マップ
  • 戦術:資産、人員、知識、計画、プロセス、プロジェクト
  • 意味:用語、定義、ルール、メタデータ、教育、ガバナンス
  • 指標:先行、遅行、兆候

テクニカルアーキテクチャー(パフォーマンスダッシュボードに直接つながるところ)

  • ディスプレイ:ダッシュボード、BIポータル、スコアカード
  • アプリケーション:モニタリング、分析、マネジメント
  • データソース:スプレッドシート、メモリーキャッシュ、DWH、データマート、レポート、ドキュメント
  • 統合:カスタムAPI、EAI(Enterprise Application Integration)、EII(Enterprise Information Integration)、クエリ実行、ETL、手動
  • データソース:レガシーシステム、パッケージのアプリ、Webページ、ファイル、サーベイ、テキスト

BIアーキテクチャー


(画像はPDFより拝借しております。)

ビジネスアーキテクチャー

・統合BI能力
 ・モニタリングレイヤー
 ・分析レイヤー
 ・レポーティングレイヤー
 ・プランニングレイヤー
・BIプラットフォーム(分析サーバ)
 ・共通のサービス、モデル、API、ファイル形式
・データデリバリーアーキテクチャー

データアーキテクチャー

Quicken Loans(アメリカの金融業者)の例

(画像はPDFより拝借しております。)

  • 企業内のソフトウェアのデータを統合し、Web経由で2日分のデータを蓄積(Real-time Store)→業務系、戦術系のダッシュボードに利用
  • Real-time Storeのデータを整形して、2ヶ月分のデータを蓄積(Operational Data Store)
  • Operational Data Storeから2週間分のデータを蓄積したものを100件ほど保持する(OLAP(online analytical processing) Cubes)→業務系、戦術系のダッシュボードに利用
  • Operational Data Storeのデータを整形して、7年分のデータを蓄積(Data Warehouse)
  • Data Warehouseのから7年分のデータを蓄積したものを250件ほど保持する(OLAP Cubes)→レポーティングや分析ツールに利用

データアーキテクチャにはいろいろあるようです。

Direct Queryアーキテクチャー

スクリーンの要素が個々のクエリに直接的にリンクしている

  • 良い点:
    • すばやくデプロイできる
    • 低コスト
  • 悪い点:
    • 浅く、ドリルダウンが制限される
    • ディメンションがない
    • ハードウェア組み込みクエリ

Query and Cacheアーキテクチャー

クエリがクエリ化可能なキャッシュとともに置かれている(In-memory or disk cache)

  • 良い点:
    • すばやくデプロイできる
    • レスポンスが速い
    • ナビゲーションが速い
  • 悪い点:
    • 静的なデータセットに縛られる

BIセマンティックレイヤー

BIツールがユーザーのためにビジネス用語で表現したクエリオブジェクトを提供

  • 良い点:
    • 抽象的なクエリオブジェクト
    • ディメンションで分けられたビュー
  • 悪い点:
    • 一般的なODBCコネクション
    • 主にDWHのヒストリカルデータ

Federated Queryアーキテクチャー

EII(Enterprise Information Integration)ツールが、スクリーンの要素と合うように複数のソースからクエリ化する

  • 良い点:
    • 複数のソース
    • セマンティックレイヤー抽出
    • デプロイが素早い
    • プロトタイプ
  • 悪い点:
    • 履歴がない
    • データの質の問題
    • 複雑性

データマートアーキテクチャー

ダッシュボードがバッチで読み込まれた永続的なデータマートに対してクエリを実行する

  • 良い点:
    • 複数のソース
    • ディメンショナルモデル
    • ヒストリカルコンテキスト
    • 素早く複雑なクエリ
  • 悪い点:
    • 即時性がない
    • 統合されていない?

Event-drivenアーキテクチャー

  • インプット:DWHや業務系システム
  • 業務系ダッシュボード:データの把握、データの集計、指標のマネジメント、イベントの検知、ルールの適用、作用/トリガー
  • アウトプット:アラート、トリガー(ワークフローエンジン)、SQL/Stored Procedures(業務系システム)

おわりに

BIダッシュボードを作成する際の洗い出しが面倒だと思っていたので、この資料で良い初期値を手に入れられました。これまでの分析業務はビジネスの一部を切り取って、疎結合なものを多く扱ってきたと思います。あるイベントの効果検証とか、ある対象の予測などです。この資料を読んで、組織の戦略などと密に絡み合い、様々な関係者の目的を成し遂げるようなBIダッシュボード作成において、組織間の調整力が強く求められるのかなと思いました。そこにデータサイエンティストの持つスキルはどうフィットするのだろうか?と思いつつ、どうやってサイエンス要素をバリューが出る形で盛り込んでやろうかと考えています。

参考情報

[1]Wayne W. Eckerson(2006). “Performance Dashboards:Measuring, Monitoring, and Managing Your Business”
[2]Sandra Durcevic(2019). “12 Best Business Intelligence Books To Get You Off the Ground With BI”, The datapine Blog

Python/Rもくもく会をプライベートで開催するための参考図書・資料をまとめる

はじめに

社内で定時後に有志で勉強会というか、その場に集まってPythonやRをもくもくと勉強をするもくもく会を開きたいと考えています。目的としては分析スキルの向上や機械学習ができるようになりたいとかいう個々人の願いを叶えることです。
色々なスキルレベルのメンバーが参加することが予想されるので、皆を幸せにするためにもレベルに応じた良い教材が必要だと思いました。
ここでは、レベルに応じて適切な教材などを忘備録として残していきたいと思います。
(私自身、全てのレベルの対象者に適切な教材を網羅しているわけではないので、随時更新していこうと思います。)

受講対象について

受講対象(PythonやRをまともに触ったことがない人)は2軸で分けるとすると以下のようになると思います。

・プログラミング経験あり/経験なし
・数学の心得あり/心得なし

  • プログラミング経験なし&数学の心得あり(アルキメデス)
    理系出身の人がメインだと思います。学部・学科によっては全然扱わないですよね。数的な思考は得意だが、それを活かすスキルが不足しているような人でしょう。眼の前におかれた数学の問題を紙とペンで解くことはできるが、仕事で使えないという感じ。私も偉そうなことは言えないですが、コードが荒れがちなので周りに良い先生がいたほうが良いと思います。

  • プログラミング経験なし&数学の心得なし(葉っぱ隊)
    一番習得に時間がかかると思います。野球やったことないのに、野球選手になりたいという人に皆さんは違和感を感じるでしょう。イメージはそんな感じです。一番時間がかかるからこそ、挫折しないための教材選びが重要かもしれません。スキル的に全裸なので、葉っぱ隊と名付けましょう。

  • プログラミング経験あり&数学の心得あり(デーサイ候補)
    最も頼もしい存在です。教科書をお渡ししておけば勝手に成長すると思います。ある程度経験を積めば分析業務を任せても良いと思います。

  • プログラミング経験あり&数学の心得なし(進捗ありマン)
    各種手法の原理を知るまではそれなりに時間がかかると思いますが、手を動かして何ができるかをすぐに味わえるので、モチベーションを維持しながら学んでいきやすいと思います。コード自体は実行できるので進捗ありマンと名付けてみましょう。

この2軸でPythonとRに関する便利な資料を探したいと思います。
ただし、どの本に関してもどのレベルの人が買っても良いとは思います。ただ、数学の心得がない中で、テイラー展開とか平均値の定理とかラグランジュ未定乗数法などの表現を目にした際に、挫折してしまう可能性があるので、適した書籍から順次広げていくのが良いと思います。なお、今回はPCでもくもくと進めれそうな書籍を選んでいます。紙とペンで進める本も重要なのですが、そのようなかた向けの書籍は取り上げていません。

アルキメデス向けの教材

Python

R

  • みんなのR 第2版
    『Rによるデータサイエンス』と迷ったのですが、プログラムの実行結果がそのまま載っている印象だったので、こちらの本がプログラミング初心者には優しいと判断しました。ほとんど数式は出てこないのですが、一般化線形モデルや時系列解析などもカバーしてくれています。また、データの前処理に関する記述もこちらの本の方が手厚いです。

葉っぱ隊向けの教材

Python

  • Pythonスタートブック [増補改訂版]
    本当にプログラミングがはじめての人向けの本です。まずはプログラミング自体に慣れたほうが良いと思います。

  • プロゲートのPython入門講座
    妻におすすめされた講座です。無料枠でもある程度学びがあるようです。環境を構築しなくても良いという点が非常に葉っぱ隊に適しているとのことです。

R

  • Rによるやさしい統計学
    Rのインストールあるいは統計学の初歩のところから、応用まで幅広く説明している本です。数式はあまり出てきませんがコードが載っているので、手を動かすことができると思います。

読み物

デーサイ候補向けの教材

Python

R

進捗ありマン向けの教材

R

  • RStudioではじめるRプログラミング入門
    プログラミング経験のある進捗ありマンであれば、R言語の扱い方をまずは知りたいだろうと思います。関数の書き方やヘルプページの使い方、オブジェクトの説明、S3の話などが詳しく書かれています。

  • 新米探偵、データ分析に挑む
    R Studioのインストール方法なども載っているので、進捗ありマンなら最初から最後まで実践できると思います。数式もほとんど出てきません。色んな分析事例をRで取り組むことで分析業務のイメージも付いてくると思います。

  • RユーザのためのRStudio[実践]入門−tidyverseによるモダンな分析フローの世界−
    R言語について何となくつかめた進捗ありマンがモダンな記法であるtidyverseを効率よく学べる良い本です。データ整形・クロス集計・可視化がモダンな記法で書けるようになると結構楽しいと思います。

Python

今後について

そもそもPythonやRに触れたことがない人にとって、Tokyo.Rの初心者セッションは少し適していないのかなと思ったので、今回は取り上げていないですが、一通り使い方をわかってもらえたら初心者セッションの資料を使ったもくもく会も開きたいと思います。最終的にはKaggle部とかを作るとかになるのかもしれませんが、そこまで行けるか行けないか。