背景
「現代数理統計学」(竹村彰通)を一通り興味ある部分はフォローしたので、実際に道具として使えるように細かい部分をフォローしていくつもりです。
今は「日本統計学会公式認定 統計検定1級対応 統計学」をぱらぱらと眺めながら章末の演習問題を解いています。Amazonでのレビューは賛否両論という感じですが、一通り別の数理統計学の教科書を読んでから入る分には簡潔にまとまっている上に、章末に理解度のチェックシートがあるので知識の抜けを補強できます。
演習問題を解いていて、確率変数の変数変換の部分の理解が甘いな、と思ったので少し真面目に手を動かしてみることにします。
問題設定
確率変数$X_0,X_1$が標準正規分布に従っているとします。このとき
Y_0 = X_0 + b X_1
のような変数変換を考えます。$f_Y(y_0)$はどのような分布に従うでしょうか。
やってみる
確率変数を$(X_0,X_1) \to (Y_0,Y_1)$に変換します。求めたいのは$f_Y(y_0)$ですので、$Y_1$の選び方には任意性があります。どうせ最後には$Y_1$ついては積分してしまいますので、計算が楽になるように選ぶ方がよいです。具体的にはヤコビアンの計算が楽になるように選びます。ここでは
Y_1 = X_1
と選びます。
x_0 = y_0-by_1,\;x_1=y_1
ですので
\begin{eqnarray}
f_X(x_0)f_X(x_1)dx_0dx_1 &=& f_X(y_0-by_1)f_X(y_1)\left|\frac{\partial (x_0,x_1)}{\partial (y_0,y_1)} \right|dy_0dy_1
\end{eqnarray}
となります。つまり、
f_Y(y_0,y_1) = f_X(y_0-by_1)f_X(y_1)\left|\frac{\partial (x_0,x_1)}{\partial (y_0,y_1)} \right|
です。面倒なのはヤコビアンの計算だけですね。今知りたいのは$f_Y(y_0)$ですので、これを$y_1$について積分してやります。
f_Y(y_0) = \int_{-\infty}^{+\infty} f_X(y_0-by_1)f_X(y_1)\left|\frac{\partial (x_0,x_1)}{\partial (y_0,y_1)} \right| dy_1 \propto e^{-\frac{y_0^2}{2(b^2+1)}}
となります。
数値実験による確認
コードは次の通りです。
# include <random>
# include <iostream>
# include <fstream>
# include <cmath>
int main()
{
std::mt19937 mt(0);
std::normal_distribution<> gauss(0, 1);
std::uniform_real_distribution<> r(0, 1);
double a = 1, b = r(mt);
int Iter = 1000000;
for (int iter = 0; iter < Iter; iter++)
{
double x0 = gauss(mt);
double x1 = gauss(mt);
double y0 = a * x0 + b * x1;
double y1 = x1;
std::cout << y0 << std::endl;
}
std::ofstream ofs("calc.txt");
for (double y0 = -10; y0 <= 10; y0 += 0.1)
{
//ofs << y0 << " " << 1. / 2. / M_PI * fabs(1. / a + b / a) * sqrt(2 * M_PI) * exp(-y0 * y0 / 2. / (b * b + 1)) / sqrt(b * b + 1)*0.1 << std::endl;;
ofs << y0 << " " << exp(-y0 * y0 / 2. / (b * b + 1)) << std::endl;;
}
return 0;
}
離散分布の場合
離散分布の場合も同じようにして計算できますが、ヤコビアンがない分わかりやすいです。きれいな形に書けるかどうかは別ですが。
条件付き確率
P(A|B) = P(A,B)/P(B)
ですので、これもあまり頭を捻らなくても同じような感じで計算できそうですね。
まとめ
これまではその都度頭を捻って計算していましたが、半機械的に計算出来そうです。最終的に積分してしまうことがわかっている確率変数はうまいこと選ぶことが重要そうです。これは教科書を斜め読みしているだけでは気が付かないポイントでした。