第四章 多次元確率変数の分布

この章で目指すこと

これまでの章では、主に1つの確率変数 $X$ の分布を見てきました。 しかし、実際の薬学・生命科学データでは、複数の量を同時に観測します。

• 血中濃度 $X$ と薬効指標 $Y$
• 遺伝子発現量 $X_1,\ldots,X_p$
• 細胞の表現型カテゴリ別の個数 $N_1,\ldots,N_k$
• 施設、患者、測定時点という階層をもつ反復データ
• 複数の潜在集団が混ざったバイオマーカー分布

多次元確率分布を学ぶと、「それぞれの変数がどう分布するか」に加えて、次を考えられるようになります。

1. 2つ以上の変数がどのように一緒に変動するか
2. ある変数の値が分かったとき、別の変数の分布がどう変わるか
3. 複数の変数を足す、割る、線形結合することで分布がどう変わるか
4. 観測されていない群や個体差を含むモデルをどう組み立てるか

この章では、初歩的な確率表から始め、統計検定1級で頻出する導出と、応用統計の医薬生物分野で重要なモデルまでつなげます。

まず「1行に複数の測定値がある」と考える

多次元という言葉は難しく見えますが、出発点は表形式データの1行です。 例えば患者1人から、次の3項目を同時に測ったとします。

患者	AUC	$C_{\max}$	奏効
A	85	6.4	1

この1行を確率変数 $(X,Y,Z)$ としてまとめて扱うのが多次元確率分布です。 用語は、表をどのように見るかに置き換えられます。

用語	表データでの直感
同時分布	複数の列を一緒に見る
周辺分布	ほかの列をいったん無視して1列だけ見る
条件付き分布	「AUCが高い患者だけ」のように条件で絞って見る
共分散・相関	2列が一緒に増減するかを要約する
変数変換	複数列から和、比、スコアなど新しい列を作る

1. 多次元確率変数

同じ標本空間 $\Omega$ 上で定義された確率変数

$$X_1,\ldots,X_p$$

をまとめて、

$$\mathbf X= \begin{pmatrix} X_1\\ \vdots\\ X_p \end{pmatrix}$$

と書き、$p$ 次元確率変数、または確率ベクトルと呼びます。

例えば、1人の患者から、

$$\mathbf X= \begin{pmatrix} \mathrm{AUC}\\ C_{\max}\\ \text{薬効スコア} \end{pmatrix}$$

を観測するなら、3次元確率変数です。 重要なのは、3つの1変量分布を別々に考えるだけでは、変数間の関係を表せないことです。

2. 同時確率分布

離散型の同時確率関数

離散型確率変数 $X,Y$ に対して、

$$p_{X,Y}(x,y)=P(X=x,Y=y)$$

を同時確率関数といいます。 次の2条件を満たします。

$$p_{X,Y}(x,y)\geq0$$ $$\sum_x\sum_y p_{X,Y}(x,y)=1$$

例えば、バイオマーカー陽性を $X=1$、陰性を $X=0$、治療反応ありを $Y=1$、なしを $Y=0$ とします。 同時確率が次の表で与えられたとします。

	$Y=0$ 反応なし	$Y=1$ 反応あり	行和
$X=0$ 陰性	0.42	0.08	0.50
$X=1$ 陽性	0.18	0.32	0.50
列和	0.60	0.40	1.00

表の各セルが同時確率です。 例えば、

$$P(X=1,Y=1)=0.32$$

は「バイオマーカー陽性かつ治療反応あり」の割合です。

連続型の同時確率密度関数

連続型確率変数 $X,Y$ に対して、関数 $f_{X,Y}(x,y)$ が、

$$f_{X,Y}(x,y)\geq0$$ $$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f_{X,Y}(x,y)\,dx\,dy=1$$

を満たし、領域 $A$ に入る確率が、

$$P\{(X,Y)\in A\} = \iint_A f_{X,Y}(x,y)\,dx\,dy$$

で与えられるとき、$f_{X,Y}$ を同時確率密度関数といいます。

1変量では曲線の下の面積が確率でした。 2変量では、密度曲面の下で、平面上の領域 $A$ の上にある体積が確率になります。

台を先に確認する

多次元分布では、式だけでなく、どの $(x,y)$ が取りうるかという台が特に重要です。 例えば、

$$f_{X,Y}(x,y)=c \qquad(0<x<y<1)$$

とします。台は単位正方形全体ではなく、直線 $y=x$ より上の三角形です。 三角形の面積は $1/2$ なので、

$$1=\int_0^1\int_0^y c\,dx\,dy =c\int_0^1y\,dy =\frac{c}{2}$$

より、

$$c=2$$

です。

同時分布関数

離散型・連続型をまとめて、

$$F_{X,Y}(x,y)=P(X\leq x,Y\leq y)$$

を同時分布関数といいます。 連続型では、

$$F_{X,Y}(x,y) =\int_{-\infty}^x\int_{-\infty}^y f_{X,Y}(u,v)\,dv\,du$$

です。十分滑らかなら、

$$f_{X,Y}(x,y) =\frac{\partial^2}{\partial x\,\partial y} F_{X,Y}(x,y)$$

として密度を復元できます。

3. 周辺分布

同時分布から一部の変数だけを見る操作を周辺化といいます。 離散型では、不要な変数について足し合わせます。

$$p_X(x)=\sum_y p_{X,Y}(x,y)$$ $$p_Y(y)=\sum_x p_{X,Y}(x,y)$$

先ほどの表では、

$$P(X=1)=0.18+0.32=0.50$$ $$P(Y=1)=0.08+0.32=0.40$$

です。

連続型では、不要な変数について積分します。

$$f_X(x)=\int_{-\infty}^{\infty}f_{X,Y}(x,y)\,dy$$ $$f_Y(y)=\int_{-\infty}^{\infty}f_{X,Y}(x,y)\,dx$$

三角形の台で $f_{X,Y}(x,y)=2$ なら、固定した $x$ に対して $y$ は $x<y<1$ を動くので、

$$f_X(x)=\int_x^1 2\,dy=2(1-x) \qquad(0<x<1)$$

です。一方、固定した $y$ に対して $x$ は $0<x<y$ を動くので、

$$f_Y(y)=\int_0^y2\,dx=2y \qquad(0<y<1)$$

です。 積分する変数だけでなく、台から積分範囲を決めることが得点の分かれ目になります。

4. 独立性

$X,Y$ が独立であるとは、すべての適切な $x,y$ について、同時分布が周辺分布の積に分解できることです。

離散型では、

$$p_{X,Y}(x,y)=p_X(x)p_Y(y)$$

連続型では、

$$f_{X,Y}(x,y)=f_X(x)f_Y(y)$$

です。 同時分布関数を使えば、

$$F_{X,Y}(x,y)=F_X(x)F_Y(y)$$

とも表せます。

先ほどのバイオマーカー表では、独立なら、

$$P(X=1,Y=1) =P(X=1)P(Y=1) =0.50\times0.40 =0.20$$

のはずです。しかし実際は0.32なので、独立ではありません。

5. 条件付き確率分布

離散型

$p_X(x)>0$ のとき、

$$p_{Y\mid X}(y\mid x) =P(Y=y\mid X=x) =\frac{p_{X,Y}(x,y)}{p_X(x)}$$

です。 バイオマーカー陽性者に限った治療反応率は、

$$P(Y=1\mid X=1) =\frac{0.32}{0.50} =0.64$$

です。全体の反応率0.40より高く、$X$ の情報を得ることで $Y$ の分布が変わっています。

連続型

$f_X(x)>0$ のとき、

$$f_{Y\mid X}(y\mid x) =\frac{f_{X,Y}(x,y)}{f_X(x)}$$

です。 三角形の例で $Y=y$ を固定すると、

$$f_{X\mid Y}(x\mid y) =\frac{2}{2y} =\frac1y \qquad(0<x<y)$$

となります。 つまり $Y=y$ が分かった後では、$X$ は区間 $(0,y)$ 上の一様分布です。

連鎖律

条件付き分布の定義を変形すると、

$$p_{X,Y}(x,y) =p_X(x)p_{Y\mid X}(y\mid x)$$

です。3変数なら、

$$p(x,y,z) =p(x)p(y\mid x)p(z\mid x,y)$$

と分解できます。 一般に、

$$p(x_1,\ldots,x_p) =p(x_1)\prod_{j=2}^p p(x_j\mid x_1,\ldots,x_{j-1})$$

です。 階層モデルやベイズモデルは、この条件付き分布の積としてモデルを組み立てます。

6. 条件付き期待値

離散型では、

$$E[Y\mid X=x] =\sum_y y\,p_{Y\mid X}(y\mid x)$$

連続型では、

$$E[Y\mid X=x] =\int_{-\infty}^{\infty} y f_{Y\mid X}(y\mid x)\,dy$$

です。

ここで大切なのは、$E[Y\mid X=x]$ は固定した $x$ に対する数ですが、

$$E[Y\mid X]$$

は $X$ の値によって変わる確率変数だという点です。

先ほどの二値データでは、

$$E[Y\mid X=0]=P(Y=1\mid X=0)=\frac{0.08}{0.50}=0.16$$ $$E[Y\mid X=1]=P(Y=1\mid X=1)=0.64$$

です。したがって、

$$E[Y\mid X] = \begin{cases} 0.16 & (X=0)\\ 0.64 & (X=1) \end{cases}$$

という確率変数になります。

反復期待値の法則

条件付き期待値をさらに平均すると、もとの期待値に戻ります。

$$\boxed{E\{E[Y\mid X]\}=E[Y]}$$

離散型で導出すると、

$$\begin{aligned} E\{E[Y\mid X]\} &=\sum_x E[Y\mid X=x]p_X(x)\\ &=\sum_x\sum_y y\,p_{Y\mid X}(y\mid x)p_X(x)\\ &=\sum_x\sum_y y\,p_{X,Y}(x,y)\\ &=\sum_y y\,p_Y(y)\\ &=E[Y] \end{aligned}$$

です。 バイオマーカー例では、

$$0.50\times0.16+0.50\times0.64=0.40$$

となり、全体の反応率と一致します。

全分散の公式

全体のばらつきは、各層の内部のばらつきと、層ごとの平均の違いに分けられます。

$$\boxed{ \mathrm{Var}(Y) =E\{\mathrm{Var}(Y\mid X)\} +\mathrm{Var}\{E[Y\mid X]\} }$$

導出では、

$$Y-E[Y] =\{Y-E[Y\mid X]\} +\{E[Y\mid X]-E[Y]\}$$

と分解して二乗します。 交差項の条件付き期待値は、

$$E\left[ \{Y-E[Y\mid X]\} \{E[Y\mid X]-E[Y]\} \mid X \right]=0$$

なので消えます。

医薬生物データでは、全体のばらつきを次のように分解できます。

全体のばらつき	条件付き分散	条件付き平均の分散
患者のAUCのばらつき	同じ遺伝子型内の個体差	遺伝子型間の平均差
ウェルの発現量のばらつき	同じプレート内の測定差	プレート間差
細胞発現量のばらつき	同じ細胞型内の差	細胞型間の平均差

この分解が、階層モデルや分散分析の基本的な考え方になります。

7. 共分散と相関係数

共分散

$X,Y$ の共分散を、

$$\mathrm{Cov}(X,Y) =E[(X-E[X])(Y-E[Y])]$$

と定義します。 展開すると、

$$\begin{aligned} \mathrm{Cov}(X,Y) &=E[XY-XE[Y]-E[X]Y+E[X]E[Y]]\\ &=E[XY]-E[X]E[Y] \end{aligned}$$

なので、

$$\boxed{\mathrm{Cov}(X,Y)=E[XY]-E[X]E[Y]}$$

です。

• 正：$X$ が大きいとき $Y$ も大きい傾向
• 負：$X$ が大きいとき $Y$ は小さい傾向
• 0：線形な連動がない

独立なら、

$$E[XY]=E[X]E[Y]$$

なので共分散は0です。 逆は一般には成り立ちません。 例えば $X$ が原点対称な分布に従い、$Y=X^2$ とすれば、強い依存関係があっても $\mathrm{Cov}(X,Y)=E[X^3]=0$ となることがあります。

相関係数

単位に依存しないよう共分散を標準化したものが相関係数です。

$$\rho_{X,Y} =\mathrm{Corr}(X,Y) =\frac{\mathrm{Cov}(X,Y)} {\sqrt{\mathrm{Var}(X)\mathrm{Var}(Y)}}$$

コーシー・シュワルツの不等式から、

$$-1\leq\rho_{X,Y}\leq1$$

です。 $|\rho|=1$ なら、確率1で $Y=aX+b$ という完全な線形関係があります。

共分散行列

$p$ 次元確率ベクトル $\mathbf X$ の平均ベクトルと共分散行列を、

$$\boldsymbol\mu=E[\mathbf X]$$ $$\boldsymbol\Sigma =E[(\mathbf X-\boldsymbol\mu) (\mathbf X-\boldsymbol\mu)^\mathsf T]$$

と定義します。 第 $(i,j)$ 成分は、

$$\Sigma_{ij}=\mathrm{Cov}(X_i,X_j)$$

です。対角成分は各変数の分散です。

任意のベクトル $\mathbf a$ に対して、

$$\mathbf a^\mathsf T\boldsymbol\Sigma\mathbf a =\mathrm{Var}(\mathbf a^\mathsf T\mathbf X) \geq0$$

共分散行列を単変量と2変量へほどく

単変量なら、$\boldsymbol\Sigma$ はただの分散 $\sigma^2$ で、$a^T\Sigma a=a^2\sigma^2=\operatorname{Var}(aX)$ です。

2変量なら、$\boldsymbol\Sigma$ の対角に $\operatorname{Var}(X_1),\operatorname{Var}(X_2)$、対角以外に $\operatorname{Cov}(X_1,X_2)$ が入ります。$\mathbf a=(a_1,a_2)^T$ とすると、二次形式は $a_1^2\operatorname{Var}(X_1)+a_2^2\operatorname{Var}(X_2)+2a_1a_2\operatorname{Cov}(X_1,X_2)$ です。

1. $\boldsymbol\Sigma\mathbf a$ を計算して各変数の分散・共分散を重み付けする。
2. 左から $\mathbf a^T$ を掛けて全成分を1個の分散へ足し合わせる。

なので、共分散行列は半正定値です。 これは共分散行列として許される数値の組み合わせに制約があることを意味します。

INTERACTIVE

多次元分布を動かして理解する

相関、和、潜在集団、乱数変換を切り替え、式が表す形を確認します。

同時・条件付き畳み込み混合分布Box–Muller

相関係数 ρ

0.60

条件 X=x

1.00

条件付き平均

条件付き分散

独立性

相関係数 $\rho$ を0に近づけると、二変量正規分布の等高線は円に近づきます。 $|\rho|$ を大きくすると楕円が細くなり、$X=x$ を固定した後の $Y$ の条件付き分散 $1-\rho^2$ は小さくなります。

8. 確率変数の和と畳み込み

離散型の畳み込み

独立な離散型確率変数 $X,Y$ の和を、

$$S=X+Y$$

とします。$S=s$ が起こる組み合わせをすべて足すと、

$$\begin{aligned} P(S=s) &=\sum_x P(X=x,Y=s-x)\\ &=\sum_x P(X=x)P(Y=s-x) \end{aligned}$$

です。したがって、

$$\boxed{ p_S(s)=\sum_x p_X(x)p_Y(s-x) }$$

となります。この演算を畳み込みといいます。

例えば、独立なベルヌーイ変数 $X,Y\sim\mathrm{Bernoulli}(p)$ の和は、

$$P(S=0)=(1-p)^2$$ $$P(S=1)=2p(1-p)$$ $$P(S=2)=p^2$$

なので、

$$S\sim\mathrm{Bin}(2,p)$$

です。

連続型の畳み込み

連続型では、$S=X+Y$ の分布関数を、

$$\begin{aligned} F_S(s) &=P(X+Y\leq s)\\ &=\int_{-\infty}^{\infty} \int_{-\infty}^{s-x} f_{X,Y}(x,y)\,dy\,dx \end{aligned}$$

と書きます。 独立なら $f_{X,Y}(x,y)=f_X(x)f_Y(y)$ なので、$s$ で微分して、

$$\boxed{ f_S(s) =\int_{-\infty}^{\infty} f_X(x)f_Y(s-x)\,dx }$$

を得ます。

母関数で和を見る

独立なら積率母関数は、

$$\begin{aligned} M_{X+Y}(t) &=E[e^{t(X+Y)}]\\ &=E[e^{tX}e^{tY}]\\ &=M_X(t)M_Y(t) \end{aligned}$$

です。この積が同じ分布族の母関数になると、再生性が分かります。

$X,Y$ の分布	必要な条件	$X+Y$ の分布
ポアソン	独立	平均母数を足したポアソン
正規	独立でなくても同時正規なら可	平均を足し、共分散を含めて分散を計算した正規
ガンマ	独立かつ尺度母数が同じ	形状母数を足したガンマ
二項	独立かつ成功確率が同じ	試行回数を足した二項

薬物曝露を複数時点の寄与の和として近似する場合や、複数区間のイベント数を足す場合に、この考え方を使います。 ただし、同一患者内の反復測定は独立でないことが多く、分散には共分散項が必要です。

$$\mathrm{Var}(X+Y) =\mathrm{Var}(X)+\mathrm{Var}(Y) +2\mathrm{Cov}(X,Y)$$

9. 多次元の変数変換とヤコビアン

なぜ行列式が必要か

1変量の変数変換では、長さの伸縮を、

$$\left|\frac{dx}{dy}\right|$$

で補正しました。 2変量では、微小な長方形の面積が変換によってどれだけ伸縮するかを、ヤコビアンの行列式で補正します。

$$\begin{pmatrix} U\\V \end{pmatrix} = \begin{pmatrix} g_1(X,Y)\\g_2(X,Y) \end{pmatrix}$$

が1対1で逆変換

$$X=x(u,v),\qquad Y=y(u,v)$$

をもつとします。このとき、

$$f_{U,V}(u,v) =f_{X,Y}\{x(u,v),y(u,v)\} \left| \frac{\partial(x,y)}{\partial(u,v)} \right|$$

です。ただし、

$$\frac{\partial(x,y)}{\partial(u,v)} = \begin{vmatrix} \dfrac{\partial x}{\partial u} & \dfrac{\partial x}{\partial v}\\ \dfrac{\partial y}{\partial u} & \dfrac{\partial y}{\partial v} \end{vmatrix}$$

です。

変数変換公式は「確率の保存」から生じる

変数変換の前後で座標の目盛りは変わっても、対応する領域に点が入る確率そのものは変わりません。

$(X,Y)$ 平面上の微小長方形の面積を $dx,dy$、対応する $(U,V)$ 平面上の微小長方形の面積を $du,dv$ とすると、

$$f_{X,Y}(x,y),dx,dy \approx f_{U,V}(u,v),du,dv$$

です。逆変換 $(x,y)=\{x(u,v),y(u,v)\}$ による面積倍率は、

$$dx,dy = \left| \frac{\partial(x,y)}{\partial(u,v)} \right| du,dv$$

なので、

$$f_{U,V}(u,v) = f_{X,Y}\{x(u,v),y(u,v)\} \left| \frac{\partial(x,y)}{\partial(u,v)} \right|$$

が得られます。

つまり、ヤコビアンは式を合わせるための飾りではなく、座標変換によって1目盛りの面積が変わる分を補い、確率を保存する係数です。

分かっている変換	密度公式で使うもの
逆変換 $(x,y)=(x(u,v),y(u,v))$	逆変換のヤコビアン行列式の絶対値
順変換 $(u,v)=(g_1(x,y),g_2(x,y))$	順変換のヤコビアン行列式の絶対値の逆数

逆数にする向きを間違えやすいため、まず「どちらの密度を求めたいか」と「どちら向きの変換を微分したか」を書きます。

和の分布を変数変換で導く

$$S=X+Y,\qquad T=X$$

とおきます。逆変換は、

$$X=T,\qquad Y=S-T$$

です。ヤコビアンは、

$$\left| \frac{\partial(x,y)}{\partial(s,t)} \right| = \left| \begin{matrix} 0&1\\ 1&-1 \end{matrix} \right| =1$$

なので、

$$f_{S,T}(s,t)=f_{X,Y}(t,s-t)$$

です。$T$ を積分して消せば、

$$f_S(s)=\int f_{X,Y}(t,s-t)\,dt$$

となり、独立なら畳み込み公式が得られます。

線形変換

$$\mathbf Y=\mathbf A\mathbf X+\mathbf b$$

なら、

$$E[\mathbf Y] =\mathbf A\boldsymbol\mu+\mathbf b$$ $$\mathrm{Cov}(\mathbf Y) =\mathbf A\boldsymbol\Sigma\mathbf A^\mathsf T$$

線形変換 $A\Sigma A^T$ の単変量対応と計算順序

単変量なら、$Y=aX+b$、$E[Y]=a\mu+b$、$\operatorname{Var}(Y)=a^2\sigma^2$ です。$\mathbf A\boldsymbol\Sigma\mathbf A^T$ は、この $a^2\sigma^2$ を複数の入力・出力へ拡張した式です。定数 $\mathbf b$ は平均を移動しますが、ばらつきは変えません。

1. $\mathbf A\mathbf X$ で元の変数を重み付き和へ組み替える。
2. $\mathbf A\boldsymbol\Sigma$ で出力側から見た共分散を作る。
3. 右から $\mathbf A^T$ を掛け、出力どうしの共分散までそろえる。

です。 この公式は、複数バイオマーカーの加重スコア、対比、主成分、線形予測量の平均と分散を計算する基本になります。

10. Box–Muller変換

Box–Muller変換は、独立な2個の一様乱数から、独立な2個の標準正規乱数を作る方法です。 単なる計算公式ではなく、次の3つを一度に確認できる重要な例です。

1. 多次元の変数変換でヤコビアンを使う理由
2. 累積分布関数の逆関数から乱数を作る方法
3. 二変量標準正規分布を極座標で見ると、半径と角度が独立になること

何を入力し、何が得られるか

$$U_1,U_2\overset{\mathrm{iid}}{\sim}U(0,1)$$

を入力とします。$\mathrm{iid}$ は、$U_1,U_2$ が同じ一様分布に従い、互いに独立であることを表します。

まず、単位正方形上の点 $(U_1,U_2)$ を、半径 $R$ と角度 $\Theta$ に変換します。

$$R=\sqrt{-2\log U_1}, \qquad \Theta=2\pi U_2$$

次に、極座標から直交座標へ戻します。

$$Z_1=R\cos\Theta, \qquad Z_2=R\sin\Theta$$

すると、

$$Z_1,Z_2\overset{\mathrm{iid}}{\sim}N(0,1)$$

となります。

入力・中間量	決めているもの	変化させると何が分かるか
$U_1$	原点からの距離 $R$	$U_1$ が0に近いほど、裾にある大きな正規乱数が生じやすい
$U_2$	方向 $\Theta$	半径を変えず、点を原点の周りに回転させる
$(Z_1,Z_2)$	正規平面上の位置	2個の独立な標準正規乱数が同時に得られる

INTERACTIVE

正方形から正規平面への変数変換

同じ点を、一様乱数、半径・角度、標準正規の3つの座標系で追跡します。

一様乱数 U₁

0.35

一様乱数 U₂

0.18

半径とヤコビアン

角度

標準正規乱数

図をどの順番で読めばよいか

図中の灰色の点は、同じ乱数標本を3つの座標系で表したものです。赤い点は、スライダーで選んだ1組の乱数です。

1. 左：$(U_1,U_2)$ は単位正方形全体に一様に分布する
2. 中：$U_1$ を半径 $R$、$U_2$ を角度 $\Theta$ へ変換する
3. 右：$(R,\Theta)$ を直交座標 $(Z_1,Z_2)$ へ戻す
4. 下：$dR\times d\Theta$ の長方形が、面積約 $R,dR,d\Theta$ の扇形へ変わる

$U_1$ だけを変えると $\Theta$ は固定されたまま赤い点が半径方向へ動きます。$U_2$ だけを変えると $R$ は固定されたまま赤い点が円周方向へ動きます。

変換ごとの役割は次の通りです。

変換	数式	座標上で起きること
一様乱数から半径	$R=\sqrt{-2\log U_1}$	正方形の横座標を原点からの距離へ変える
一様乱数から角度	$\Theta=2\pi U_2$	正方形の縦座標を一周の方向へ変える
極座標から直交座標	$Z_1=R\cos\Theta$	半径の横方向成分を取り出す
極座標から直交座標	$Z_2=R\sin\Theta$	半径の縦方向成分を取り出す

下段の扇形は、ヤコビアン、

$$\left| \frac{\partial(z_1,z_2)} {\partial(r,\theta)} \right| =r$$

の幾何学的な意味です。同じ $dr,d\theta$ でも、原点から遠いほど扇形の面積が大きくなるため、倍率 $r$ が必要になります。

導出1：目標となる二変量標準正規分布

独立な標準正規確率変数 $Z_1,Z_2$ の1変量密度は、

$$\phi(z)=\frac{1}{\sqrt{2\pi}}e^{-z^2/2}$$

です。独立なら同時密度は積になるため、

$$\begin{aligned} f_{Z_1,Z_2}(z_1,z_2) &=\phi(z_1)\phi(z_2)\\ &=\frac{1}{2\pi} \exp\left\{-\frac{z_1^2+z_2^2}{2}\right\} \end{aligned}$$

です。指数部分が $z_1^2+z_2^2$ だけで決まるため、密度は原点からの距離が等しい円周上で同じ値になります。 これが極座標を使う理由です。

導出2：直交座標を極座標へ変換する

極座標変換を、

$$z_1=r\cos\theta,\qquad z_2=r\sin\theta$$

とします。台は、

$$r>0,\qquad 0\leq\theta<2\pi$$

です。

高校数学の1変数変換では $|dx/dy|$ で長さを補正しました。2変数では、微小面積の倍率をヤコビアンで補正します。

$$\frac{\partial(z_1,z_2)}{\partial(r,\theta)} = \begin{vmatrix} \dfrac{\partial z_1}{\partial r} & \dfrac{\partial z_1}{\partial\theta}\\[4pt] \dfrac{\partial z_2}{\partial r} & \dfrac{\partial z_2}{\partial\theta} \end{vmatrix}$$

各偏微分を計算すると、

$$\frac{\partial z_1}{\partial r}=\cos\theta, \qquad \frac{\partial z_1}{\partial\theta}=-r\sin\theta$$ $$\frac{\partial z_2}{\partial r}=\sin\theta, \qquad \frac{\partial z_2}{\partial\theta}=r\cos\theta$$

なので、

$$\begin{aligned} \frac{\partial(z_1,z_2)}{\partial(r,\theta)} &= \begin{vmatrix} \cos\theta&-r\sin\theta\\ \sin\theta&r\cos\theta \end{vmatrix}\\ &=r\cos^2\theta+r\sin^2\theta\\ &=r \end{aligned}$$

です。したがってヤコビアンの絶対値も $r$ です。

この $r$ は、半径が大きくなるほど同じ幅 $d\theta$ の扇形の面積が大きくなることを表します。実際、微小面積は、

$$dz_1\,dz_2=r\,dr\,d\theta$$

と変換されます。

導出3：半径と角度の同時密度を求める

多次元の変数変換公式より、

$$\begin{aligned} f_{R,\Theta}(r,\theta) &=f_{Z_1,Z_2}(r\cos\theta,r\sin\theta) \left| \frac{\partial(z_1,z_2)}{\partial(r,\theta)} \right|\\ &=\frac{1}{2\pi} \exp\left[-\frac{r^2(\cos^2\theta+\sin^2\theta)}{2}\right]r\\ &=\frac{1}{2\pi}re^{-r^2/2} \end{aligned}$$

です。ただし、$r>0,\ 0\leq\theta<2\pi$ です。

ここで、

$$f_{R,\Theta}(r,\theta) =\underbrace{re^{-r^2/2}}_{f_R(r)} \underbrace{\frac{1}{2\pi}}_{f_\Theta(\theta)}$$

と、$r$ だけの関数と $\theta$ だけの関数の積に分かれます。したがって、

$$R\ \text{と}\ \Theta\ \text{は独立}$$

です。

周辺密度を積分から確認すると、

$$\begin{aligned} f_R(r) &=\int_0^{2\pi}f_{R,\Theta}(r,\theta)\,d\theta\\ &=\int_0^{2\pi}\frac{1}{2\pi}re^{-r^2/2}\,d\theta\\ &=re^{-r^2/2} \end{aligned}$$

であり、

$$\begin{aligned} f_\Theta(\theta) &=\int_0^\infty f_{R,\Theta}(r,\theta)\,dr\\ &=\frac{1}{2\pi} \int_0^\infty re^{-r^2/2}\,dr\\ &=\frac{1}{2\pi} \end{aligned}$$

です。最後の積分では $w=r^2/2$ とおくと $dw=r\,dr$ となり、

$$\int_0^\infty re^{-r^2/2}\,dr =\int_0^\infty e^{-w}\,dw =1$$

を使いました。

導出4：半径を一様乱数から作る

$R$ の累積分布関数を、密度から計算します。$r>0$ に対して、

$$\begin{aligned} F_R(r) &=P(R\leq r)\\ &=\int_0^r se^{-s^2/2}\,ds\\ &=\left[-e^{-s^2/2}\right]_0^r\\ &=1-e^{-r^2/2} \end{aligned}$$

です。

逆関数法では、一様乱数 $V\sim U(0,1)$ に対して、

$$V=F_R(R)=1-e^{-R^2/2}$$

とおいて $R$ について解きます。

$$e^{-R^2/2}=1-V$$

両辺の自然対数を取ると、

$$-\frac{R^2}{2}=\log(1-V)$$

したがって、$R>0$ を使って、

$$R=\sqrt{-2\log(1-V)}$$

です。$V$ が一様分布なら $1-V$ も一様分布なので、$1-V$ を改めて $U_1$ と書けば、

$$R=\sqrt{-2\log U_1}$$

を得ます。

実装では $U_1=0$ だと $\log0$ が定義できないため、必ず、

$$0<U_1<1$$

となる一様乱数を使います。連続一様分布では端点を取る確率は0ですが、コンピュータ上の乱数では明示的な確認が必要です。

導出5：角度を一様乱数から作る

$\Theta$ は区間 $[0,2\pi)$ 上の一様分布なので、累積分布関数は、

$$F_\Theta(\theta)=\frac{\theta}{2\pi} \qquad(0\leq\theta<2\pi)$$

です。$U_2=F_\Theta(\Theta)$ とおけば、

$$\Theta=2\pi U_2$$

となります。

$U_1$ と $U_2$ は独立なので、それぞれから作られる $R$ と $\Theta$ も独立です。この $(R,\Theta)$ を直交座標へ戻したものが、

$$\boxed{ \begin{aligned} Z_1&=\sqrt{-2\log U_1}\cos(2\pi U_2)\\ Z_2&=\sqrt{-2\log U_1}\sin(2\pi U_2) \end{aligned} }$$

です。

なぜ $Z_1,Z_2$ は「独立な」標準正規分布になるのか

ここまでで作った $(R,\Theta)$ の同時密度は、二変量標準正規分布を極座標へ変換したときの同時密度と一致しています。 したがって、直交座標へ戻した $(Z_1,Z_2)$ の同時密度は、

$$f_{Z_1,Z_2}(z_1,z_2) =\frac{1}{2\pi}e^{-(z_1^2+z_2^2)/2}$$

です。さらに、

$$\frac{1}{2\pi}e^{-(z_1^2+z_2^2)/2} = \left(\frac{1}{\sqrt{2\pi}}e^{-z_1^2/2}\right) \left(\frac{1}{\sqrt{2\pi}}e^{-z_2^2/2}\right)$$

と周辺密度の積に分解できます。よって、

$$Z_1\sim N(0,1),\qquad Z_2\sim N(0,1),\qquad Z_1\perp Z_2$$

です。

単に $\mathrm{Cov}(Z_1,Z_2)=0$ だから独立なのではありません。同時密度が周辺密度の積に分解できることが独立性の根拠です。

発展：変換全体のヤコビアンで直接確かめる

Box–Muller変換を、

$$(U_1,U_2) \longrightarrow (R,\Theta) \longrightarrow (Z_1,Z_2)$$

という2つの変換の合成として見ます。

まず、

$$R=\sqrt{-2\log U_1}$$

の両辺を2乗すると、

$$R^2=-2\log U_1$$

です。$U_1$ で微分すると、

$$2R\frac{dR}{dU_1} =-\frac{2}{U_1}$$

なので、

$$\frac{dR}{dU_1} =-\frac{1}{RU_1}$$

です。また、

$$\frac{d\Theta}{dU_2}=2\pi$$

です。$R$ は $U_2$ に依存せず、$\Theta$ は $U_1$ に依存しないため、

$$\begin{aligned} \left| \frac{\partial(r,\theta)} {\partial(u_1,u_2)} \right| &= \left| \begin{matrix} -1/(ru_1)&0\\ 0&2\pi \end{matrix} \right|\\ &=\frac{2\pi}{ru_1} \end{aligned}$$

です。

次の極座標変換のヤコビアンは、

$$\left| \frac{\partial(z_1,z_2)} {\partial(r,\theta)} \right| =r$$

でした。合成変換のヤコビアンは積になるため、

$$\begin{aligned} \left| \frac{\partial(z_1,z_2)} {\partial(u_1,u_2)} \right| &= \left| \frac{\partial(z_1,z_2)} {\partial(r,\theta)} \right| \left| \frac{\partial(r,\theta)} {\partial(u_1,u_2)} \right|\\ &=r\times\frac{2\pi}{ru_1}\\ &=\frac{2\pi}{u_1} \end{aligned}$$

です。

求めたいのは $(Z_1,Z_2)$ の密度なので、逆向きのヤコビアンを使います。

$$\left| \frac{\partial(u_1,u_2)} {\partial(z_1,z_2)} \right| =\frac{u_1}{2\pi}$$

入力 $(U_1,U_2)$ の同時密度は単位正方形上で1です。また、

$$u_1 =e^{-r^2/2} =e^{-(z_1^2+z_2^2)/2}$$

なので、

$$\begin{aligned} f_{Z_1,Z_2}(z_1,z_2) &= f_{U_1,U_2}(u_1,u_2) \left| \frac{\partial(u_1,u_2)} {\partial(z_1,z_2)} \right|\\ &=1\times\frac{u_1}{2\pi}\\ &=\frac{1}{2\pi} e^{-(z_1^2+z_2^2)/2} \end{aligned}$$

となり、二変量標準正規密度が直接得られます。

この導出では、極座標の $r$ と $dR/dU_1$ に含まれる $1/r$ が打ち消し合う点が重要です。

数値例：$U_1=e^{-1},\ U_2=1/8$

半径と角度は、

$$R=\sqrt{-2\log(e^{-1})}=\sqrt2, \qquad \Theta=2\pi\times\frac18=\frac{\pi}{4}$$

です。したがって、

$$Z_1=\sqrt2\cos\frac{\pi}{4}=1, \qquad Z_2=\sqrt2\sin\frac{\pi}{4}=1$$

となります。また、

$$Z_1^2+Z_2^2=1^2+1^2=2=R^2$$

であり、極座標の関係も確認できます。

発展：半径の二乗とカイ二乗分布

$$R^2=Z_1^2+Z_2^2$$

です。独立な標準正規乱数の平方和なので、

$$R^2\sim\chi^2_2$$

です。一方、Box–Mullerの式から直接、

$$R^2=-2\log U_1$$

とも書けます。$S=R^2$ とすると、

$$P(S\leq s)=1-e^{-s/2}$$

なので、$S$ は率 $1/2$ の指数分布でもあります。これは、

$$\chi^2_2=\mathrm{Gamma}(1,\text{scale }2)$$

という分布の関係と一致します。

アルゴリズムとして書く

Pythonで式をそのまま実装すると、次のようになります。

from math import cos, log, pi, sin, sqrt


def box_muller(u1: float, u2: float) -> tuple[float, float]:
    if not 0.0 < u1 < 1.0:
        raise ValueError("u1 must satisfy 0 < u1 < 1")
    if not 0.0 <= u2 < 1.0:
        raise ValueError("u2 must satisfy 0 <= u2 < 1")

    radius = sqrt(-2.0 * log(u1))
    angle = 2.0 * pi * u2
    z1 = radius * cos(angle)
    z2 = radius * sin(angle)
    return z1, z2

ここで log は常用対数ではなく自然対数です。1組の $(U_1,U_2)$ から2個の正規乱数が得られるため、$Z_1$ だけ使って $Z_2$ を捨てる必要はありません。

発展：一般の正規分布・相関した正規分布へ広げる

標準正規乱数 $Z$ から平均 $\mu$、標準偏差 $\sigma$ の正規乱数を作るには、

$$X=\mu+\sigma Z$$

とします。

多変量正規乱数では、Box–Muller変換などで独立な標準正規ベクトル $\mathbf Z$ を作り、共分散行列を、

$$\boldsymbol\Sigma=\mathbf L\mathbf L^\mathsf T$$

とCholesky分解して、

$$\mathbf X=\boldsymbol\mu+\mathbf L\mathbf Z$$

と変換します。すると、

$$E[\mathbf X]=\boldsymbol\mu, \qquad \mathrm{Cov}(\mathbf X)=\boldsymbol\Sigma$$

Cholesky変換を単変量から読む

単変量なら、$\sigma^2=\sigma\cdot\sigma$ と分解し、$X=\mu+\sigma Z$ とするだけです。多変量の $\mathbf L$ は標準偏差 $\sigma$ の役割を持ちますが、各成分を伸縮するだけでなく、混ぜ合わせて相関も作ります。

1. 独立標準正規ベクトル $\mathbf Z$ を作る。
2. 下三角行列 $\mathbf L$ を掛けて分散と相関を入れる。
3. 平均ベクトル $\boldsymbol\mu$ を足して中心を移す。

です。

医薬・生命科学での使い方

使用場面	正規乱数が表すもの	Box–Mullerから先に必要なモデル
母集団薬物動態シミュレーション	個体間変動、残差変動	変動の分散・共分散、分布仮定、共変量モデル
検定力・標本サイズのシミュレーション	仮想患者の測定誤差や反応値	効果量、脱落、群間差、解析手法
バイオマーカー予測	相関した複数測定値	共分散行列、測定誤差、外部検証条件
モンテカルロ積分	確率モデルからの標本	推定対象、重み、収束診断

例えば、母集団薬物動態モデルの個体間変動を、

$$\boldsymbol\eta\sim N_p(\mathbf0,\boldsymbol\Omega)$$

と置く場合、独立標準正規乱数を $\boldsymbol\Omega$ に対応する線形変換へ通して $\boldsymbol\eta$ を生成できます。 ただし、正規乱数を生成できることと、実データに正規分布を仮定してよいことは別問題です。分布仮定、外れ値、個体内・個体間階層、推定した共分散行列の妥当性を確認する必要があります。

よくある間違い

1. $U_1=0$ を許す：$\log0$ が定義できません。
2. 常用対数を使う：式の $\log$ は自然対数です。
3. ヤコビアンの $r$ を落とす：面積補正がなくなり、$f_R(r)$ を正しく導けません。
4. $U_1,U_2$ の独立性を書かない：$R$ と $\Theta$ の独立性を保証する前提です。
5. 台を書かない：$r>0,\ 0\leq\theta<2\pi$ が必要です。
6. 無相関だけで独立と結論する：この場合の独立性は、同時密度の積への分解から示します。

Box–Muller変換は、モンテカルロ法、検定統計量のシミュレーション、薬物動態モデルの仮想患者生成などで、正規乱数を作る原理の1つです。 実務の乱数ライブラリでは別の高速アルゴリズムが使われることもありますが、変数変換、逆関数法、ヤコビアンをまとめて理解する教材として非常に重要です。

11. 多項分布

二項分布から多項分布へ

1回の試行の結果が $k$ 個のカテゴリのいずれかに入り、カテゴリ $j$ の確率を $p_j$ とします。

$$p_j\geq0, \qquad \sum_{j=1}^k p_j=1$$

$n$ 回の独立な試行で、カテゴリ $j$ に入った回数を $N_j$ とすると、

$$N_1+\cdots+N_k=n$$

です。この個数ベクトルが、

$$(N_1,\ldots,N_k) \sim\mathrm{Multinomial}(n;p_1,\ldots,p_k)$$

に従うといいます。

確率関数は、

$$P(N_1=n_1,\ldots,N_k=n_k) = \frac{n!}{n_1!\cdots n_k!} \prod_{j=1}^k p_j^{n_j}$$

です。ただし、

$$n_j=0,1,\ldots,n, \qquad \sum_{j=1}^k n_j=n$$

です。

多項係数が現れる理由

特定の並び、例えば最初の $n_1$ 回がカテゴリ1、次の $n_2$ 回がカテゴリ2という1つの並びの確率は、

$$\prod_{j=1}^k p_j^{n_j}$$

です。 同じ個数 $(n_1,\ldots,n_k)$ を生む並び方は、

$$\frac{n!}{n_1!\cdots n_k!}$$

通りあるため、両者を掛けます。

多項定理

多項定理は、

$$(x_1+\cdots+x_k)^n = \sum_{\substack{n_1+\cdots+n_k=n\\n_j\geq0}} \frac{n!}{n_1!\cdots n_k!} \prod_{j=1}^k x_j^{n_j}$$

です。 $x_j=p_j$ とおけば、

$$(p_1+\cdots+p_k)^n=1$$

なので、多項分布の確率を全個数パターンについて足すと1になることが分かります。

平均・分散・共分散の導出

試行 $i$ がカテゴリ $j$ なら1、それ以外なら0となる指示変数を、

$$I_{ij} = \begin{cases} 1 & (\text{試行 }i\text{ がカテゴリ }j)\\ 0 & (\text{それ以外}) \end{cases}$$

とします。すると、

$$N_j=\sum_{i=1}^n I_{ij}$$

です。各 $I_{ij}$ はベルヌーイ分布なので、

$$E[N_j]=np_j$$ $$\mathrm{Var}(N_j)=np_j(1-p_j)$$

です。

$j\neq\ell$ について、同じ試行が2つのカテゴリに同時に入ることはないので、

$$I_{ij}I_{i\ell}=0$$

です。したがって、

$$\begin{aligned} \mathrm{Cov}(I_{ij},I_{i\ell}) &=E[I_{ij}I_{i\ell}] -E[I_{ij}]E[I_{i\ell}]\\ &=-p_jp_\ell \end{aligned}$$

となり、独立な試行間の共分散は0なので、

$$\boxed{ \mathrm{Cov}(N_j,N_\ell) =-np_jp_\ell \qquad(j\neq\ell) }$$

です。 1カテゴリの個数が増えれば、合計 $n$ が固定されているため、別カテゴリの個数は減ります。この制約が負の共分散を生みます。

周辺化・統合・条件付け

1カテゴリだけ見れば、

$$N_j\sim\mathrm{Bin}(n,p_j)$$

です。 複数カテゴリをまとめた個数も、確率を足した二項分布になります。 例えば、

$$N_A+N_B \sim\mathrm{Bin}(n,p_A+p_B)$$

です。

また、カテゴリ1と2の合計が $m$ だと分かったとき、

$$N_1\mid(N_1+N_2=m) \sim \mathrm{Bin}\left( m,\frac{p_1}{p_1+p_2} \right)$$

です。

医薬生物分野では、細胞を複数表現型に分類した個数、患者の有害事象グレード、遺伝子型の個数、シーケンスリードのカテゴリ別個数などに現れます。 ただし、細胞やリードに過分散・クラスタリングがある場合、単純な多項分布では分散を過小評価することがあります。

12. 多変量正規分布

$p$ 次元確率ベクトル $\mathbf X$ が平均ベクトル $\boldsymbol\mu$、正定値共分散行列 $\boldsymbol\Sigma$ をもつ多変量正規分布に従うとき、

$$\mathbf X\sim N_p(\boldsymbol\mu,\boldsymbol\Sigma)$$

と書きます。 密度は、

$$f(\mathbf x) = \frac{1} {(2\pi)^{p/2}|\boldsymbol\Sigma|^{1/2}} \exp\left[ -\frac12 (\mathbf x-\boldsymbol\mu)^\mathsf T \boldsymbol\Sigma^{-1} (\mathbf x-\boldsymbol\mu) \right]$$

多変量正規密度の行列を1変量へ戻す

$p=1$ なら、$|\boldsymbol\Sigma|^{1/2}=\sigma$、$\boldsymbol\Sigma^{-1}=1/\sigma^2$ なので、通常の正規密度 $\exp[-(x-\mu)^2/(2\sigma^2)]/(\sqrt{2\pi}\sigma)$ になります。

$\mathbf x-\boldsymbol\mu$ は「平均との差」、$\boldsymbol\Sigma^{-1}$ は「分散と相関による標準化」、左の転置ベクトルは「各方向の寄与を足して1個の数へ戻す」操作です。行列式 $|\boldsymbol\Sigma|$ は分布の楕円全体の広がりを表します。

です。

マハラノビス距離

何を測る距離なのか

平均ベクトル $\boldsymbol\mu$ から観測点 $\mathbf x$ までの差を、

$$\boldsymbol\delta=\mathbf x-\boldsymbol\mu$$

とします。マハラノビス距離の二乗は、

$$D^2 \equiv (\mathbf x-\boldsymbol\mu)^\mathsf T \boldsymbol\Sigma^{-1} (\mathbf x-\boldsymbol\mu) =\boldsymbol\delta^\mathsf T \boldsymbol\Sigma^{-1}\boldsymbol\delta$$

で定義されます。距離そのものは、

$$D=\sqrt{D^2}$$

マハラノビス距離の掛け算を順番に追う

単変量なら、$D^2=(x-\mu)^2/\sigma^2$、すなわち標準得点の二乗です。

1. $\boldsymbol\delta=\mathbf x-\boldsymbol\mu$ で平均との差を作る。
2. $\boldsymbol\Sigma^{-1}\boldsymbol\delta$ で、ばらつきが大きい方向を弱く、相関に逆らう方向を強く補正する。
3. $\boldsymbol\delta^T$ を左から掛けて、補正後の各成分を1個の距離の二乗へ合計する。
4. 距離が必要なら最後に平方根を取る。

です。

通常のユークリッド距離、

$$d_E(\mathbf x,\boldsymbol\mu) =\sqrt{ (\mathbf x-\boldsymbol\mu)^\mathsf T (\mathbf x-\boldsymbol\mu) }$$

は、すべての方向を同じ尺度で測ります。一方、マハラノビス距離は、次の2点を補正します。

1. 分散の違い：ばらつきの大きい変数の1単位と、ばらつきの小さい変数の1単位を同じ重さにしない
2. 相関：よく一緒に動く方向は珍しくないと評価し、相関に逆らう方向は珍しいと評価する

したがって、マハラノビス距離は「物理的に何単位離れたか」ではなく、その分布にとって何標準偏差ぶん不自然な方向へ離れたかを測る距離です。

相関がない場合：標準化したユークリッド距離

変数が無相関で、

$$\boldsymbol\Sigma = \begin{pmatrix} \sigma_1^2&0\\ 0&\sigma_2^2 \end{pmatrix}$$

なら、逆行列は、

$$\boldsymbol\Sigma^{-1} = \begin{pmatrix} 1/\sigma_1^2&0\\ 0&1/\sigma_2^2 \end{pmatrix}$$

です。したがって、

$$D^2 = \frac{(x_1-\mu_1)^2}{\sigma_1^2} + \frac{(x_2-\mu_2)^2}{\sigma_2^2}$$

となります。

つまり、各変数を、

$$z_j=\frac{x_j-\mu_j}{\sigma_j}$$

と標準化した後のユークリッド距離です。分散が100の変数での1単位差は小さく、分散が0.01の変数での1単位差は非常に大きい、と区別できます。

2変量で相関を入れた式を導く

標準偏差を $\sigma_1,\sigma_2$、相関係数を $\rho$ とすると、共分散行列は、

$$\boldsymbol\Sigma = \begin{pmatrix} \sigma_1^2&\rho\sigma_1\sigma_2\\ \rho\sigma_1\sigma_2&\sigma_2^2 \end{pmatrix}$$

です。行列式は、

$$|\boldsymbol\Sigma| =\sigma_1^2\sigma_2^2(1-\rho^2)$$

なので、$|\rho|<1$ なら逆行列が存在します。2次正方行列の逆行列公式から、

$$\boldsymbol\Sigma^{-1} = \frac{1} {\sigma_1^2\sigma_2^2(1-\rho^2)} \begin{pmatrix} \sigma_2^2&-\rho\sigma_1\sigma_2\\ -\rho\sigma_1\sigma_2&\sigma_1^2 \end{pmatrix}$$

です。

$$a=\frac{x_1-\mu_1}{\sigma_1}, \qquad b=\frac{x_2-\mu_2}{\sigma_2}$$

と標準化すると、

$$\boxed{ D^2 = \frac{a^2-2\rho ab+b^2}{1-\rho^2} }$$

となります。

中央の、

$$-2\rho ab$$

が相関補正です。$\rho>0$ のとき、$a,b$ が同符号なら距離を小さくし、異符号なら距離を大きくします。

白色化すると普通の距離になる

相関を除き、各方向の分散を1へ揃える変換を白色化といいます。2変量では、例えば、

$$w_1=a$$ $$w_2 =\frac{b-\rho a}{\sqrt{1-\rho^2}}$$

と変換できます。このとき、

$$\begin{aligned} w_1^2+w_2^2 &=a^2+\frac{(b-\rho a)^2}{1-\rho^2}\\ &=\frac{a^2-2\rho ab+b^2}{1-\rho^2}\\ &=D^2 \end{aligned}$$

です。

つまり、原空間のマハラノビス距離は、白色化空間でのユークリッド距離です。

INTERACTIVE

ユークリッド距離とマハラノビス距離

分散・相関・観測点を動かし、同じ点が白色化後の空間でどこへ移るかを比較します。

標準偏差 σ₁

1.40

標準偏差 σ₂

0.80

相関 ρ

0.75

観測点 x₁

1.80

観測点 x₂

-0.80

ユークリッド距離

マハラノビス距離 D

χ²₂累積確率

図の原空間では等距離線が楕円ですが、白色化後は円になります。赤い点を動かしたとき、右側の原点から赤い点までの距離がマハラノビス距離 $D$ です。

面白い現象：ユークリッド距離が同じでも珍しさは全く違う

平均が $(0,0)$、標準偏差がともに1、相関係数が $\rho=0.8$ の二変量正規分布を考えます。

$$\mathbf x_A=(2,2)^\mathsf T, \qquad \mathbf x_B=(2,-2)^\mathsf T$$

とすると、どちらもユークリッド距離は、

$$d_E(\mathbf x_A,\mathbf0) =d_E(\mathbf x_B,\mathbf0) =\sqrt{2^2+2^2} =\sqrt8$$

です。

しかし、マハラノビス距離は異なります。$\mathbf x_A$ では、

$$\begin{aligned} D_A^2 &=\frac{2^2-2(0.8)(2)(2)+2^2}{1-0.8^2}\\ &=\frac{1.6}{0.36}\\ &\approx4.44 \end{aligned}$$

です。一方、$\mathbf x_B$ では、

$$\begin{aligned} D_B^2 &=\frac{2^2-2(0.8)(2)(-2)+(-2)^2}{1-0.8^2}\\ &=\frac{14.4}{0.36}\\ &=40 \end{aligned}$$

です。

正の相関が強い集団では「両方高い」は相関方向に沿った変化ですが、「一方が高く他方が低い」は相関に逆らう珍しい組み合わせです。 各変数を別々に見ると、どちらも平均から2標準偏差しか離れていません。それでも、$\mathbf x_B$ は周辺的には極端でなくても、組み合わせとして非常に異常です。

これはマハラノビス距離の最も面白く、実用的な性質です。

確率楕円との関係

$$D^2=c$$

を満たす点は、2次元では楕円、3次元では楕円体を作ります。 共分散行列を固有値分解して、

$$\boldsymbol\Sigma =\mathbf Q\boldsymbol\Lambda\mathbf Q^\mathsf T$$

と書くと、$\mathbf Q$ の列である固有ベクトルが楕円の向き、固有値 $\lambda_j$ が各方向の分散を表します。

$$D^2=c$$

の楕円の半軸長は、

$$\sqrt{c\lambda_1}, \qquad \sqrt{c\lambda_2}$$

固有値分解を単変量と図形で読む

単変量なら、共分散は $\sigma^2$ という1個の数で、固有方向は数直線そのもの、固有値は $\sigma^2$ です。多変量では、$\mathbf Q^T$ で楕円の軸に合わせて座標を回転し、$\boldsymbol\Lambda$ で各軸方向の分散を読み、$\mathbf Q$ で元の向きへ戻します。

計算上は $\mathbf Q\boldsymbol\Lambda\mathbf Q^T$ を右から読みます。点を固有軸へ移す、軸ごとに伸縮する、元の座標へ戻す、という3段階です。

です。分散が大きい固有方向では楕円が長くなり、同じ座標差でも距離への寄与は小さくなります。

多変量正規分布で、母平均と母共分散行列が既知なら、

$$D^2\sim\chi^2_p$$

です。したがって、確率 $\gamma$ を含む楕円領域は、

$$D^2\leq\chi^2_{p,\gamma}$$

で作れます。$p=2$ では、

含まれる確率 $\gamma$	$\chi^2_{2,\gamma}$	楕円の境界
0.50	1.386	$D^2=1.386$
0.95	5.991	$D^2=5.991$
0.99	9.210	$D^2=9.210$

です。インタラクティブ図の破線は50%領域、実線は95%領域を表します。

正規密度における「驚き度」

多変量正規分布の対数密度を取ると、

$$\begin{aligned} -2\log f(\mathbf x) &=p\log(2\pi) +\log|\boldsymbol\Sigma| +D^2 \end{aligned}$$

です。平均と共分散行列を固定すれば、最初の2項は定数なので、$D^2$ が大きい点ほど密度が低くなります。

この意味でマハラノビス距離は、分布モデルから見た驚き度として解釈できます。ただし、$D^2$ が大きいことだけで、測定ミス、病態変化、新しいサブタイプのどれかまでは判定できません。

2点間の距離と単位変換への不変性

平均からの距離だけでなく、2点 $\mathbf x,\mathbf y$ の距離を、

$$D(\mathbf x,\mathbf y) = \sqrt{ (\mathbf x-\mathbf y)^\mathsf T \boldsymbol\Sigma^{-1} (\mathbf x-\mathbf y) }$$

と定義できます。

また、可逆な線形変換、

$$\mathbf z=\mathbf A\mathbf x+\mathbf b$$

を行い、共分散行列も、

$$\boldsymbol\Sigma_z =\mathbf A\boldsymbol\Sigma\mathbf A^\mathsf T$$

と変換すれば、マハラノビス距離は変わりません。 mg/Lをng/mLへ変えるような単位換算や、座標軸の回転によって異常度が変わらないのが利点です。

医薬・生命科学で何が分かるか

観測情報	マハラノビス距離で学べること	注意点
AUCと$C_{\max}$	各指標は正常範囲でも、組み合わせが典型的PK関係から外れていないか	用量、体重、投与経路などで層別が必要
複数バイオマーカー	単独値では見逃す多変量外れ値	病型や治療群を混ぜた共分散は不適切なことがある
フローサイトメトリーの複数マーカー	既知細胞集団からの距離、QC異常	非正規性、補償、ゲーティング、バッチ差
scRNA-seqの低次元表現	参照集団から離れた細胞状態	細胞を独立反復とみなさず、個体・バッチを考慮
分析法・製造QCの複数特性	個別規格内でも全体として不自然なロット	管理限界は事前に定義し、性能評価する
化合物記述子	既知化学空間からの外挿度	距離が活性や安全性を直接保証するわけではない

Topic：マハラノビス距離と毒性予測

最初に区別する2つの問い

QSARや機械学習による毒性予測では、次の2つは別の問いです。

1. 毒性予測：その化合物が毒性陽性となる確率や毒性値はいくつか
2. 適用領域：その化合物は、モデルが学習した化学空間の内側にあるか

毒性モデルが出す、

$$\widehat P(\text{毒性}\mid\mathbf x)$$

は1つ目への答えです。一方、マハラノビス距離、

$$D^2 =(\mathbf x-\bar{\mathbf x}_{\mathrm{train}})^{\mathsf T} \mathbf S_{\mathrm{train}}^{-1} (\mathbf x-\bar{\mathbf x}_{\mathrm{train}})$$

は2つ目への答えとして使えます。

ここで、

• $\mathbf x$：候補化合物の記述子ベクトル
• $\bar{\mathbf x}_{\mathrm{train}}$：学習化合物の記述子平均
• $\mathbf S_{\mathrm{train}}$：学習化合物の記述子共分散行列

です。

したがって、

$$\boxed{ \text{マハラノビス距離は毒性の強さではなく、学習化学空間からの外れ方を測る} }$$

と理解します。

INTERACTIVE TOPIC

毒性確率と適用領域を分けて見る

仮想QSARの出力と、候補化合物が学習済み化学空間の内側にあるかを同時に確認します。

候補化合物 x₁

1.10

候補化合物 x₂

0.40

仮想モデル P（毒性）

マハラノビス距離 D²

適用領域の判定

図の紫破線は、説明用に設定した仮想毒性モデルの判別境界です。青い楕円は、記述子が二変量正規分布に近いと仮定した95%化学空間です。

橙色の候補点を動かすと、毒性確率と $D^2$ は別々に変化します。

• 毒性確率が高く、95%領域内：モデルは学習領域内で毒性を予測している
• 毒性確率が高く、95%領域外：毒性懸念は無視できないが、数値は外挿なので追加検証が必要
• 毒性確率が低く、95%領域内：学習領域内で低毒性と予測している
• 毒性確率が低く、95%領域外：低い予測値を安全性の根拠にしてはいけない

図の毒性確率は説明用であり、実在する化合物やエンドポイントを予測した値ではありません。

何の情報から何を学べるか

入力情報	マハラノビス距離で分かること	それだけでは分からないこと
分子量、logP、極性表面積など	学習済み物性空間からの外挿度	毒性機序、代謝活性化
構造記述子・埋め込み表現	学習化合物との多変量的な近さ	特定部分構造の毒性寄与
複数in vitro assay応答	既知応答パターンからの外れ方	ヒトでの用量反応や曝露
トランスクリプトーム署名	参照毒性応答からの距離	臓器毒性の因果機序
毒性予測モデルの学習データ	補間か外挿かの目安	外部検証性能、確率校正

距離は「似ているか」を多変量的に要約しますが、「なぜ毒性が生じるか」や「ヒトで安全か」を直接示しません。

毒性予測での基本手順

1. 毒性エンドポイントを明確にする
2. 学習データだけで記述子の標準化、平均、共分散行列を求める
3. 学習データで毒性予測モデルを構築する
4. 交差検証・外部検証で予測性能と確率校正を評価する
5. 候補化合物について毒性予測値とマハラノビス距離を別々に計算する
6. 適用領域外なら、類似化合物、構造アラート、代謝、追加assayを確認する

平均や共分散を全データから計算してから訓練・評価データへ分割すると、評価データの情報が学習へ漏れるデータリーケージになります。適用領域の計算も学習データだけで定義します。

数値例：高い毒性確率でも外挿の場合

標準化した2つの記述子が相関 $\rho=0.65$ をもつとします。候補化合物を、

$$\mathbf x=(2,-2)^{\mathsf T}$$

とすると、

$$\begin{aligned} D^2 &=\frac{2^2-2(0.65)(2)(-2)+(-2)^2} {1-0.65^2}\\ &=\frac{4+5.2+4}{0.5775}\\ &=\frac{13.2}{0.5775}\\ &\approx22.86 \end{aligned}$$

です。

二変量正規分布の95%境界 $5.991$ と比べると、

$$22.86>5.991$$

なので、この候補は95%化学空間の外側です。

図の仮想モデル、

$$\operatorname{logit} \{P(\text{毒性})\} =-0.3+0.9x_1-0.6x_2$$

では、

$$\begin{aligned} \eta &=-0.3+0.9(2)-0.6(-2)\\ &=2.7 \end{aligned}$$

なので、

$$P(\text{毒性}) =\frac{1}{1+e^{-2.7}} \approx0.937$$

です。

毒性確率は約93.7%ですが、外挿予測なので「93.7%という数値を高精度に信頼できる」とは結論できません。一方、安全側の判断では、この結果を単に捨てるのでもなく、確認試験や類似化合物調査を優先する警告として扱います。

発展：線形QSARのレバレッジとの関係

$n$ 個の学習化合物について、中心化した $n\times p$ 記述子行列を $\mathbf X_c$ とします。標本共分散行列は、

$$\mathbf S =\frac{1}{n-1}\mathbf X_c^{\mathsf T}\mathbf X_c$$

です。

中心化した候補記述子を $\mathbf x_c$ とすると、

$$\begin{aligned} D^2 &=\mathbf x_c^{\mathsf T} \mathbf S^{-1} \mathbf x_c\\ &=(n-1) \mathbf x_c^{\mathsf T} (\mathbf X_c^{\mathsf T}\mathbf X_c)^{-1} \mathbf x_c \end{aligned}$$

です。切片を含む線形モデルでの候補点のレバレッジは、

$$h =\frac1n +\mathbf x_c^{\mathsf T} (\mathbf X_c^{\mathsf T}\mathbf X_c)^{-1} \mathbf x_c$$

なので、

$$\boxed{ h=\frac1n+\frac{D^2}{n-1} }$$

となります。

したがって、学習重心から遠い化合物ほどレバレッジも大きくなります。Williams plotなどで使われる、

$$h^*=\frac{3(p+1)}{n}$$

は警告レバレッジの経験的基準ですが、普遍的な安全境界ではありません。

毒性予測で特に注意すること

1. 距離が小さくても予測が正しいとは限らない：モデルのバイアス、ラベル誤り、クラス不均衡、確率校正不良は別問題です。
2. 距離が大きくても毒性とは限らない：新規な化学構造や物性を示すだけで、毒性の方向は示しません。
3. 楕円内にデータの空白があり得る：多峰性・非線形な化学空間では、重心距離だけでは局所的なデータ密度を表せません。
4. フィンガープリントにはそのまま適用しにくい：連続記述子の共分散に基づく距離なので、構造類似度や近傍法も併用します。
5. $p\geq n$ では逆行列が不安定：記述子選択、PCA、縮小共分散、正則化を検討します。
6. 代謝物・反応性・曝露を忘れない：親化合物の記述子空間だけでは、代謝活性化や臓器曝露を十分に表せないことがあります。
7. 適用領域は目的別に定義する：肝毒性、心毒性、変異原性など、エンドポイントが変われば学習空間と信頼性も変わります。

OECDの(Q)SAR評価では、モデルの適用領域と限界を明確にし、領域外の予測を外挿として扱うことが重視されています。ただし、信頼できる領域と信頼できない領域の間に絶対的な一本の境界があるわけではありません。

参考：

• OECD (Q)SAR Assessment Framework
• OECD Guidance on the validation of (Q)SAR models
• NCATS Tox21 Program Goals

PKで起こる「各項目は正常だが組み合わせが変」

AUCと$C_{\max}$ が正に相関する参照集団を考えます。

• AUCも$C_{\max}$も高い患者は、相関方向に沿うため距離がそれほど大きくない場合がある
• AUCが高いのに$C_{\max}$が低い患者は、各値が単独では基準内でも、参照関係から外れて大きな距離になる場合がある

後者は、吸収速度、採血時点、服薬状況、測定誤差、製剤差などを確認するきっかけになります。ただし、距離だけで原因を断定してはいけません。

集団を混ぜると異常判定が逆転することがある

健康群と疾患群、異なる細胞型、異なる投与量群を1つの楕円でまとめると、群間差が共分散を大きくします。 その結果、本当の外れ値が楕円内へ隠れるマスキングや、正常な少数群が外れ値扱いされるスワンピングが起こります。

「どの参照集団からの距離か」を先に定義することが重要です。

発展：Hotellingの $T^2$ と判別分析

標本平均 $\bar{\mathbf X}$ が仮説平均 $\boldsymbol\mu_0$ からどれだけ離れているかを、

$$T^2 =n(\bar{\mathbf X}-\boldsymbol\mu_0)^\mathsf T \mathbf S^{-1} (\bar{\mathbf X}-\boldsymbol\mu_0)$$

Hotellingの $T^2$ は単変量の何に当たるか

単変量なら、$T^2=n(\bar X-\mu_0)^2/S^2=\{\sqrt n(\bar X-\mu_0)/S\}^2$ で、1標本 $t$ 統計量の二乗です。多変量では、平均差を標本共分散の逆行列で標準化し、相関を考慮した平均差の大きさへ変えています。

で測るのがHotellingの $T^2$ 統計量です。1標本t検定を多変量へ拡張した形で、マハラノビス距離に標本数 $n$ を掛けています。

また、群ごとの平均ベクトルは異なるが共分散行列は共通と仮定すると、「各群平均へのマハラノビス距離が最も小さい群を選ぶ」という考え方が線形判別分析につながります。

実務での重要な注意点

1. 共分散行列を同じデータから推定した場合：$D^2\sim\chi^2_p$ をそのまま使えないことがあります。有限標本補正やHotellingの $T^2$ を検討します。
2. 外れ値が共分散を膨らませる場合：外れ値自身が距離を小さく見せるマスキングが起こります。ロバスト共分散推定を検討します。
3. $p\geq n$ または強い多重共線性：$\mathbf S$ が特異になり逆行列を計算できません。変数選択、PCA、縮小推定、正則化を使います。
4. 高次元で距離が集中する場合：変数数が多いほど、近い点と遠い点の差が相対的に小さくなります。次元削減と外部検証が必要です。
5. 非正規・多峰性のデータ：カイ二乗閾値と楕円解釈が合わない場合があります。変換、混合モデル、ロバスト法、経験分布を検討します。
6. 距離と因果を混同しない：大きな距離は「参照分布と合わない」ことを示しますが、その原因は示しません。

線形結合は正規分布

任意の定数ベクトル $\mathbf a$ に対して、

$$\mathbf a^\mathsf T\mathbf X \sim N\left( \mathbf a^\mathsf T\boldsymbol\mu, \mathbf a^\mathsf T\boldsymbol\Sigma\mathbf a \right)$$

です。 さらに行列 $\mathbf A$ とベクトル $\mathbf b$ に対して、

$$\mathbf A\mathbf X+\mathbf b \sim N_q\left( \mathbf A\boldsymbol\mu+\mathbf b, \mathbf A\boldsymbol\Sigma\mathbf A^\mathsf T \right)$$

ベクトル $a$ と行列 $A$ の役割の違い

$\mathbf a^T\mathbf X$ は複数変数を1個の値へまとめます。単変量なら $aX$ です。一方、$\mathbf A\mathbf X$ は複数の重み付き和を同時に作り、別のベクトルへ変換します。行列 $\mathbf A$ の各行が、1個ずつ新しい変数を作る重みです。

分散は、1出力なら $\mathbf a^T\boldsymbol\Sigma\mathbf a$、複数出力なら $\mathbf A\boldsymbol\Sigma\mathbf A^T$ と覚えると、式の形を取り違えにくくなります。

です。

周辺分布

多変量正規分布の一部を取り出した周辺分布も正規分布です。 平均ベクトルと共分散行列を、

$$\boldsymbol\mu= \begin{pmatrix} \boldsymbol\mu_1\\ \boldsymbol\mu_2 \end{pmatrix}, \qquad \boldsymbol\Sigma= \begin{pmatrix} \boldsymbol\Sigma_{11}&\boldsymbol\Sigma_{12}\\ \boldsymbol\Sigma_{21}&\boldsymbol\Sigma_{22} \end{pmatrix}$$

と分割すれば、

$$\mathbf X_1\sim N(\boldsymbol\mu_1,\boldsymbol\Sigma_{11})$$

です。

条件付き分布

多変量正規分布では条件付き分布も正規分布です。

$$\mathbf X_1\mid\mathbf X_2=\mathbf x_2 \sim N(\boldsymbol\mu_{1\mid2},\boldsymbol\Sigma_{1\mid2})$$

ただし、

$$\boldsymbol\mu_{1\mid2} = \boldsymbol\mu_1 +\boldsymbol\Sigma_{12} \boldsymbol\Sigma_{22}^{-1} (\mathbf x_2-\boldsymbol\mu_2)$$ $$\boldsymbol\Sigma_{1\mid2} = \boldsymbol\Sigma_{11} -\boldsymbol\Sigma_{12} \boldsymbol\Sigma_{22}^{-1} \boldsymbol\Sigma_{21}$$

条件付き正規分布を2変量のスカラー式へ戻す

$X_1=X$、$X_2=Y$ の2変量で、標準偏差を $\sigma_X,\sigma_Y$、相関を $\rho$ とすると、条件付き平均は $\mu_X+\rho(\sigma_X/\sigma_Y)(y-\mu_Y)$、条件付き分散は $\sigma_X^2(1-\rho^2)$ です。

1. $y-\mu_Y$ で観測された $Y$ のずれを求める。
2. $\Sigma_{22}^{-1}$ で $Y$ 自身のばらつきを標準化する。
3. $\Sigma_{12}$ で $X$ と $Y$ の連動分だけ $X$ の平均を修正する。
4. 分散から説明できた部分 $\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}$ を引く。

です。

標準化された2変量正規分布、

$$\begin{pmatrix}X\\Y\end{pmatrix} \sim N_2\left[ \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}1&\rho\\\rho&1\end{pmatrix} \right]$$

なら、

$$Y\mid X=x \sim N(\rho x,1-\rho^2)$$

です。 相関が強いほど $X$ の情報から $Y$ を狭い範囲に予測できます。

医薬生物分野での使い方

• AUC、$C_{\max}$、クリアランスの対数値を同時にモデル化する
• 複数バイオマーカーの同時分布と予測区間を作る
• 線形混合モデルのランダム効果に多変量正規分布を置く
• 遺伝子発現の低次元表現や判別分析の基礎分布として使う

現実のオミクスデータは歪み、外れ値、群の混在を含むため、多変量正規性を仮定する前に散布図、変換、ロバスト性を確認します。

13. 混合分布

まず「1個の観測値」が生まれる場面を考える

混合分布は、性質の異なる複数の集団が混ざっているが、各観測がどの集団に由来するかは直接見えないときのモデルです。

例えば、薬剤処理後の100個の細胞について応答量 $X$ を測定したとします。

• 薬剤感受性細胞では応答量が低い分布になる
• 薬剤耐性細胞では応答量が高い分布になる
• しかし、測定した応答量だけから各細胞の表現型を完全には判別できない

このとき、1個の細胞の値は次の順序で生じると考えます。

1. まず、その細胞が感受性群か耐性群かを確率的に選ぶ
2. 選ばれた群の分布から応答量を1つ生成する
3. 観測者には応答量 $X$ は見えるが、本当の群 $Z$ は見えない

ここで $Z$ を潜在クラスと呼びます。「潜在」は、存在を仮定するが直接は観測していない、という意味です。

潜在クラスを使って数式にする

潜在クラスを $Z\in\{1,\ldots,K\}$ とし、

$$P(Z=k)=\pi_k, \qquad \sum_{k=1}^K\pi_k=1$$

とします。$Z=k$ のときの条件付き密度を $f_k(x)$ とすれば、全確率の公式から、

$$\begin{aligned} f_X(x) &=\sum_{k=1}^K P(Z=k)f_{X\mid Z}(x\mid k)\\ &=\sum_{k=1}^K\pi_k f_k(x) \end{aligned}$$

です。これを有限混合分布といいます。

2つの正規分布を混ぜる場合は、

$$X\sim \pi N(\mu_1,\sigma_1^2) +(1-\pi)N(\mu_2,\sigma_2^2)$$

と表します。これは略記であり、正確には、

$$P(Z=1)=\pi, \qquad P(Z=2)=1-\pi$$ $$X\mid Z=1\sim N(\mu_1,\sigma_1^2)$$ $$X\mid Z=2\sim N(\mu_2,\sigma_2^2)$$

という二段階のモデルです。

INTERACTIVE

混合分布と階層モデルの違い

見えない集団の混在と、既知のグループに入れ子になった観測を、生成過程から比較します。

混合分布階層モデル

集団1の割合 π

0.60

集団間の距離 Δ

2.80

観測値 x

0.50

P（集団1｜x）

P（集団2｜x）

混合平均

「混合分布」タブでは、破線が各成分の重み付き密度、

$$\pi f_1(x), \qquad (1-\pi)f_2(x)$$

で、青い実線がその和 $f_X(x)$ です。観測値 $x$ を動かすと、その値がどちらの集団から生じた可能性が高いかも変わります。

「階層モデル」タブでは、全体平均から群平均、群平均から観測値が生じる順序を示します。灰色点が観測値、橙色のひし形が生の群別平均、緑色点が全体情報を借りた部分プーリング後の推定値です。詳しい導出は第14節で扱います。

観測値から潜在クラスを逆算する

生成時には「$Z$ を選んでから $X$ を生成」します。一方、解析時には $X=x$ を観測した後で、

$$P(Z=k\mid X=x)$$

を求めたいことがあります。ベイズの公式から、

$$\boxed{ P(Z=k\mid X=x) = \frac{\pi_k f_k(x)} {\sum_{j=1}^K\pi_j f_j(x)} }$$

です。この確率を、混合モデルでは所属確率や責任度と呼ぶことがあります。

分子は「成分 $k$ が、その観測値 $x$ を生じさせる重み付き密度」、分母は「すべての成分を合わせた密度」です。

例えば、

$$P(Z=1)=0.7, \qquad P(Z=2)=0.3$$ $$X\mid Z=1\sim N(0,1), \qquad X\mid Z=2\sim N(4,1)$$

とします。$x=3$ を観測したとき、

$$\begin{aligned} P(Z=1\mid X=3) &= \frac{0.7\phi(3)} {0.7\phi(3)+0.3\phi(-1)}\\ &\approx \frac{0.7(0.00443)} {0.7(0.00443)+0.3(0.24197)}\\ &\approx0.041 \end{aligned}$$

です。したがって、$x=3$ は成分2から生じた可能性が高いと推定されます。

ただし、所属確率が高くてもクラスを実験的に同定したことにはなりません。表面マーカー、遺伝子型、機能試験などによる独立した検証が必要です。

混合分布の平均を導く

成分 $k$ の平均を、

$$\mu_k=E[X\mid Z=k]$$

とします。反復期待値の法則を使うと、

$$\begin{aligned} E[X] &=E\{E[X\mid Z]\}\\ &=\sum_{k=1}^K P(Z=k)E[X\mid Z=k]\\ &=\sum_{k=1}^K\pi_k\mu_k \end{aligned}$$

です。全体平均は、成分平均の重み付き平均になります。

2成分なら、

$$\mu =\pi\mu_1+(1-\pi)\mu_2$$

です。

混合分布の分散を導く

全分散の公式、

$$\mathrm{Var}(X) =E\{\mathrm{Var}(X\mid Z)\} +\mathrm{Var}\{E[X\mid Z]\}$$

を使います。

第1項は、

$$E\{\mathrm{Var}(X\mid Z)\} =\sum_{k=1}^K\pi_k\sigma_k^2$$

であり、各成分内部のばらつきの平均です。

第2項は、

$$\mathrm{Var}\{E[X\mid Z]\} =\sum_{k=1}^K\pi_k(\mu_k-\mu)^2$$

であり、成分平均どうしの違いです。したがって、

$$\boxed{ \mathrm{Var}(X) = \underbrace{ \sum_{k=1}^K\pi_k\sigma_k^2 }_{\text{成分内分散}} + \underbrace{ \sum_{k=1}^K\pi_k(\mu_k-\mu)^2 }_{\text{成分間分散}} }$$

となります。

数値例：全体の分散が大きい理由を分ける

確率0.7で $N(0,1)$、確率0.3で $N(4,1)$ から値が生成されるとします。

全体平均は、

$$\mu=0.7(0)+0.3(4)=1.2$$

です。成分内分散は、

$$0.7(1)+0.3(1)=1$$

です。成分間分散は、

$$\begin{aligned} &0.7(0-1.2)^2+0.3(4-1.2)^2\\ &=0.7(1.44)+0.3(7.84)\\ &=1.008+2.352\\ &=3.36 \end{aligned}$$

です。したがって、

$$\mathrm{Var}(X)=1+3.36=4.36$$

です。各成分の分散は1しかありません。全体の大きな分散4.36の主因は、成分内部が不安定だからではなく、平均の異なる集団をまとめたことです。

山の数と成分数は同じとは限らない

2成分混合分布でも、成分平均が近い、片方の割合が小さい、各成分の分散が大きい、といった場合には山が1つに見えます。

逆に、1つの連続分布でも歪みや外れ値によって複数の山のように見えることがあります。したがって、

$$\text{ヒストグラムの山が2つ} \quad\not\Rightarrow\quad \text{真の生物学的集団が2つ}$$

です。

医薬・生命科学での混合分布

観測値	想定する潜在成分	学べること	追加検証
薬効スコア	レスポンダー・非レスポンダー	成分比、群別薬効分布	バイオマーカーによる再現
細胞の薬剤応答	感受性・耐性表現型	耐性細胞の割合	クローン化、機能試験
遺伝子発現量	複数細胞型・細胞状態	潜在状態の分布	マーカー、空間情報
PK指標	代謝型・非代謝型	曝露分布の異質性	遺伝型、併用薬、肝腎機能
有害事象発生時刻	早期・遅延期群	異なる発生機序の候補	臨床背景、追跡設計

例えば受容体活性評価で最大反応 $E_{\max}$ が二峰性に見えても、すぐに「受容体高発現群と低発現群」と断定してはいけません。細胞数、トランスフェクション効率、プレート、測定日、シグナル飽和など、連続的・技術的要因でも同様の形が生じます。

14. 階層モデル

まず「誰の中で測ったか」を残す

階層モデルは、観測値が個体、施設、バッチ、ウェルなどのグループに入れ子になっているときに使います。

例えば4匹のマウスから、それぞれ10枚の切片を作り、各切片で蛍光強度を測ったとします。

$$4\text{匹}\times10\text{切片}=40\text{観測}$$

ですが、独立なマウスが40匹いるわけではありません。同じマウスの切片は、遺伝背景、投与量、組織状態などを共有するため似やすくなります。

添字を、

$$Y_{ij} = \text{マウス }i\text{ の切片 }j\text{ の測定値}$$

と書くと、階層が見えます。

• $i$：独立に近い実験単位であるマウス
• $j$：同じマウス内の反復観測

階層モデルの重要な問いは、「40個の値があるか」ではなく、どの観測が同じ上位単位を共有しているかです。

生成過程を3段階で書く

最も基本的な正規階層モデルは、

$$\mu \longrightarrow \theta_i \longrightarrow Y_{ij}$$

という順序で考えます。

1. 集団全体の平均 $\mu$ がある
2. 個体ごとの平均 $\theta_i$ が集団平均の周囲にばらつく
3. 各個体の反復測定 $Y_{ij}$ が個体平均の周囲にばらつく

数式では、

$$\theta_i\mid\mu,\tau^2 \sim N(\mu,\tau^2)$$ $$Y_{ij}\mid\theta_i,\sigma^2 \sim N(\theta_i,\sigma^2)$$

と書きます。

• $\tau^2$：個体間・施設間・バッチ間の分散
• $\sigma^2$：同じ個体・施設・バッチ内の分散

別の書き方では、

$$Y_{ij}=\mu+u_i+\varepsilon_{ij}$$ $$u_i\sim N(0,\tau^2), \qquad \varepsilon_{ij}\sim N(0,\sigma^2)$$

です。$u_i$ は群 $i$ に共通するずれ、$\varepsilon_{ij}$ は各観測固有のずれです。

「階層モデル」タブでは、灰色点が観測値、橙色のひし形が群別平均、緑色点が部分プーリング後の推定値です。

同じ群の観測が似る理由

期待値は、

$$\begin{aligned} E[Y_{ij}] &=E[\mu+u_i+\varepsilon_{ij}]\\ &=\mu \end{aligned}$$

です。分散は、$u_i$ と $\varepsilon_{ij}$ が独立なら、

$$\begin{aligned} \mathrm{Var}(Y_{ij}) &=\mathrm{Var}(u_i)+\mathrm{Var}(\varepsilon_{ij})\\ &=\tau^2+\sigma^2 \end{aligned}$$

です。

同じ群 $i$ の異なる2観測 $Y_{ij},Y_{ik}$ の共分散を計算します。

$$Y_{ij}=\mu+u_i+\varepsilon_{ij}$$ $$Y_{ik}=\mu+u_i+\varepsilon_{ik}$$

両方に共通する確率変数は $u_i$ だけなので、

$$\begin{aligned} \mathrm{Cov}(Y_{ij},Y_{ik}) &=\mathrm{Cov}(u_i+\varepsilon_{ij},u_i+\varepsilon_{ik})\\ &=\mathrm{Var}(u_i)\\ &=\tau^2 \end{aligned}$$

です。

したがって、同じ群内の2観測の相関係数は、

$$\boxed{ \mathrm{ICC} = \frac{\tau^2}{\tau^2+\sigma^2} }$$

です。ICCは全分散のうち、上位グループの違いで説明される割合です。

例えば、

$$\tau=0.6, \qquad \sigma=0.8$$

なら、

$$\mathrm{ICC} =\frac{0.6^2}{0.6^2+0.8^2} =\frac{0.36}{1.00} =0.36$$

です。全体のばらつきの36%がマウス間差、64%がマウス内差に対応します。

異なる群 $i\neq h$ なら $u_i$ と $u_h$ を共有しないため、モデル上は、

$$\mathrm{Cov}(Y_{ij},Y_{hk})=0$$

です。

疑似反復がなぜ危険か

ICCが正なのに40切片を40個の独立標本として解析すると、同じマウス由来切片の重複情報を別々の情報として数えてしまいます。

各群の標本数が $m$ でICCが $\rho$ のとき、単純化したデザイン効果は、

$$\mathrm{DE}=1+(m-1)\rho$$

です。

例えば1匹当たり10切片、ICCが0.36なら、

$$\mathrm{DE} =1+(10-1)(0.36) =4.24$$

です。見かけ上40観測あっても、同程度の情報量をもつ独立標本数の目安は、

$$n_{\mathrm{eff}} \approx\frac{40}{4.24} \approx9.43$$

です。これは近似ですが、「切片数を増やすこと」と「マウス数を増やすこと」が同じではないと分かります。

3つの推定法：完全・非・部分プーリング

複数施設の反応率や複数患者のPKパラメータを推定するとき、3つの考え方があります。

方法	推定の考え方	長所	弱点
完全プーリング	全群を同じとみなし、全体平均だけを使う	安定する	本当の群間差を消す
非プーリング	各群を完全に別々に推定する	群差を保つ	小標本群が極端になりやすい
部分プーリング	群別平均と全体平均を情報量で重み付けする	群差と安定性を両立	階層分布の仮定が必要

階層モデルの代表的な利点が部分プーリングです。

部分プーリングの式を導く

群 $i$ に $n_i$ 個の観測があり、その標本平均を $\bar Y_i$ とします。

$$Y_{ij}\mid\theta_i \sim N(\theta_i,\sigma^2)$$

なら、標本平均の分布は、

$$\bar Y_i\mid\theta_i \sim N\left(\theta_i,\frac{\sigma^2}{n_i}\right)$$

です。また、群平均の集団分布を、

$$\theta_i\sim N(\mu,\tau^2)$$

とします。

$\mu,\tau^2,\sigma^2$ を既知とし、正規分布どうしを組み合わせると、$\theta_i$ の条件付き平均は、

$$E[\theta_i\mid\bar Y_i] = \frac{ \dfrac{n_i}{\sigma^2}\bar Y_i + \dfrac{1}{\tau^2}\mu }{ \dfrac{n_i}{\sigma^2} + \dfrac{1}{\tau^2} }$$

です。分子・分母に $\sigma^2\tau^2$ を掛けると、

$$\begin{aligned} E[\theta_i\mid\bar Y_i] &= \frac{ n_i\tau^2\bar Y_i+\sigma^2\mu }{ n_i\tau^2+\sigma^2 }\\ &= B_i\bar Y_i+(1-B_i)\mu \end{aligned}$$

となります。ただし、

$$\boxed{ B_i = \frac{n_i\tau^2} {n_i\tau^2+\sigma^2} }$$

です。

• $n_i$ が大きい：$B_i$ が1に近づき、群自身の平均を重視する
• $\tau^2$ が大きい：本当に群差が大きいので、群自身の平均を重視する
• $\sigma^2$ が大きい：群内データが不安定なので、全体平均を重視する

数値例：小標本施設は強く縮む

全体平均 $\mu=2$、施設平均 $\bar Y_i=5$、群間分散 $\tau^2=1$、群内分散 $\sigma^2=4$ とします。

$n_i=3$ なら、

$$B_i =\frac{3(1)}{3(1)+4} =\frac37$$

なので、

$$\begin{aligned} E[\theta_i\mid\bar Y_i] &=\frac37(5)+\frac47(2)\\ &=\frac{23}{7}\\ &\approx3.29 \end{aligned}$$

です。生の施設平均5は、全体平均2の方向へ3.29まで縮みます。

一方、$n_i=30$ なら、

$$B_i =\frac{30}{30+4} \approx0.882$$

なので、

$$E[\theta_i\mid\bar Y_i] \approx0.882(5)+0.118(2) \approx4.65$$

です。標本数が多いため、施設自身の平均を強く残します。

混合分布との関係

一般的な階層モデルを、

$$Y_i\mid\theta_i\sim p(y_i\mid\theta_i)$$ $$\theta_i\sim p(\theta_i\mid\eta)$$

とします。 $\theta_i$ を積分して見えなくすると、

$$p(y_i\mid\eta) =\int p(y_i\mid\theta_i) p(\theta_i\mid\eta)\,d\theta_i$$

となり、周辺分布は混合分布です。 つまり、混合分布と階層モデルは数学的に切り離されたものではありません。条件付き分布の階層を、潜在変数を消して外側から見ると混合分布になることがあります。

正規階層モデルなら、

$$\theta_i\sim N(\mu,\tau^2)$$ $$Y_{ij}\mid\theta_i\sim N(\theta_i,\sigma^2)$$

より、1観測だけの周辺分布は、

$$Y_{ij}\sim N(\mu,\tau^2+\sigma^2)$$

です。しかし、周辺分布だけを見ると「どの観測が同じ個体に属するか」という情報が消えます。階層モデルでは、その所属情報を残して群内相関を扱います。

混合分布と階層モデルを見分ける

問い	混合分布	階層モデル
各観測の所属は分かるか	通常は分からない	患者ID、施設IDなどが分かる
見えない変数	成分ラベル $Z$ が典型	群ごとの効果 $\theta_i$ が典型
主な目的	潜在集団の割合と成分分布を知る	群内変動と群間変動を分ける
データの形	出所不明の値が混ざる	観測が群の中に入れ子になる
代表例	レスポンダーと非レスポンダー	患者内の反復採血
典型的な誤用	山をそのまま生物学的群と断定	細胞や切片を独立標本として扱う

次の問いで判断すると分かりやすくなります。

1. 各観測の患者ID・マウスID・施設IDは既知か
2. 同じ上位単位から複数回測定しているか
3. 未知のサブタイプを推定したいのか
4. 個体差・施設差を分散として分離したいのか

未知のサブタイプが中心なら混合分布、既知の入れ子構造が中心なら階層モデルが第一候補です。両方が同時に存在するモデルもあります。

ベータ二項モデル

施設ごとに有害事象発生確率が異なる状況を考えます。

$$P_i\sim\mathrm{Beta}(\alpha,\beta)$$ $$Y_i\mid P_i=p_i\sim\mathrm{Bin}(n_i,p_i)$$

とします。

• 上位階層：施設ごとの発生確率 $P_i$ が施設間でばらつく
• 下位階層：施設 $i$ では、その施設固有の確率 $P_i$ で事象が生じる

$P_i$ を積分して消すと、$Y_i$ はベータ二項分布に従います。

$$P(Y_i=y) = \binom {n_i}y \frac{B(\alpha+y,\beta+n_i-y)} {B(\alpha,\beta)}$$

平均を $\mu=\alpha/(\alpha+\beta)$、集中度を $\kappa=\alpha+\beta$ とすると、

$$E[Y_i]=n_i\mu$$ $$\mathrm{Var}(Y_i) =n_i\mu(1-\mu) \times\frac{\kappa+n_i}{\kappa+1}$$

です。 通常の二項分布の分散 $n_i\mu(1-\mu)$ に、

$$\frac{\kappa+n_i}{\kappa+1}$$

という1以上の係数が掛かります。施設ごとの確率が同じでないため、単純な二項分布よりばらつきが大きくなります。これが過分散です。

同じ群に属する2つのベルヌーイ結果の相関は、

$$\rho=\frac{1}{\kappa+1}$$

です。$\kappa$ が小さいほど群間の不均一性が大きく、同じ群内の観測が似やすくなります。

ポアソン・ガンマ混合

$$Y_i\mid\Lambda_i=\lambda_i \sim\mathrm{Poisson}(\lambda_i)$$ $$\Lambda_i \sim\mathrm{Gamma}(\alpha,\text{rate } \beta)$$

とします。培養ウェルごとにコロニー発生率 $\Lambda_i$ が異なり、その発生率のもとでコロニー数 $Y_i$ が生じる、という二段階です。

$\Lambda_i$ を積分すると、

$$P(Y_i=y) = \frac{\Gamma(\alpha+y)} {\Gamma(\alpha)y!} \left(\frac{\beta}{\beta+1}\right)^\alpha \left(\frac{1}{\beta+1}\right)^y$$

となり、負の二項分布です。

平均は反復期待値から、

$$\begin{aligned} E[Y_i] &=E\{E[Y_i\mid\Lambda_i]\}\\ &=E[\Lambda_i]\\ &=\frac{\alpha}{\beta} \end{aligned}$$

です。分散は全分散の公式から、

$$\begin{aligned} \mathrm{Var}(Y_i) &=E\{\mathrm{Var}(Y_i\mid\Lambda_i)\} +\mathrm{Var}\{E[Y_i\mid\Lambda_i]\}\\ &=E[\Lambda_i]+\mathrm{Var}(\Lambda_i)\\ &=\frac{\alpha}{\beta} +\frac{\alpha}{\beta^2} \end{aligned}$$

なので、分散が平均を上回ります。 第1項がポアソン分布そのもののばらつき、第2項がウェル・個体・細胞間の発生率の違いです。

医薬・生命科学での階層

上位単位	下位観測	分けたい変動	典型的な解析
患者	血中濃度の時系列	個体間PK差・測定誤差	母集団PK、非線形混合効果
施設	患者の転帰	施設差・患者差	一般化線形混合モデル
マウス	切片・視野・細胞	個体差・画像内変動	線形混合モデル
ドナー	単一細胞	ドナー差・細胞差	pseudobulk、階層モデル
プレート	ウェル	プレート差・ウェル差	ランダム効果、バッチ補正
実験日	GPCR応答曲線	日間差・測定誤差	非線形階層モデル

15. 医薬生物分野で何を学べるか

観測する情報	モデル化する分布	学べること	主な注意点
バイオマーカーと薬効	同時・条件付き分布	マーカー値を知った後の反応分布	交絡、選択バイアス
AUCと$C_{\max}$	多変量正規・対数変換	曝露指標の共変動、同時予測	外れ値、非線形性
表現型別細胞数	多項分布	各表現型割合とカテゴリ間共分散	過分散、細胞の独立性
患者別反復測定	階層モデル	個体内変動と個体間変動	欠測、時間相関
レスポンダー混在	混合分布	潜在群の割合と群別分布	成分数、検証可能性
コロニー・発現カウント	ポアソン・ガンマ階層	平均発生率と不均一性	ゼロ過剰、ライブラリサイズ
シミュレーション	Box–Muller・多変量正規	仮想データ、検定力、予測分布	乱数種、モデル仮定

scRNA-seqでの注意

細胞を独立な標本として扱うと、同一個体由来細胞の相関を無視します。 個体が実験単位なら、個体内に細胞が入れ子になった階層を意識する必要があります。 細胞型割合は合計1の組成データなので、カテゴリ間に機械的な負の依存が生じます。 多項分布・ディリクレ多項分布・ロジスティック正規分布などが候補になります。

薬物動態での注意

同一患者の濃度時系列は、同じクリアランスや分布容積を共有するため相関します。 母集団薬物動態では、個体パラメータを集団分布から生じるランダム効果として扱い、測定誤差と個体間変動を分けます。 単純な相関係数だけでなく、時間、用量、共変量、個体階層を含むモデルが必要です。

情報科学・機械学習との接続

機械学習の入力は、1標本につき1個の数値ではなく、通常は多数の特徴量を並べたベクトルです。

$$\boldsymbol X=(X_1,\ldots,X_p)^T$$

したがって、共分散行列、条件付き分布、線形変換、混合分布、階層構造は、そのまま多変量データ解析の基礎になります。

多変量統計の概念	情報科学・機械学習での利用	医薬・生命科学の例
共分散行列	特徴量の重複情報を把握	遺伝子群、分子記述子、PK指標
固有ベクトル	PCAによる次元削減	scRNA-seqの低次元表現
マハラノビス距離	多変量異常検知、適用領域	毒性予測の分布外化合物
条件付き分布	欠測補完、条件付き予測	バイオマーカーから薬効を予測
混合分布	ソフトクラスタリング	レスポンダー・非レスポンダー混在
階層モデル	群・施設・個体差の分離	患者内反復、バッチ、ドナー差
多変量変数変換	潜在表現、生成モデル	分子・細胞の埋め込み表現

PCAは共分散行列の固有値問題である

中心化したベクトル $\boldsymbol X$ を、長さ1の方向 $\boldsymbol a$ へ射影すると、

$$Z=\boldsymbol a^T\boldsymbol X$$

です。その分散は、

$$\operatorname{Var}(Z) =\boldsymbol a^T\boldsymbol\Sigma\boldsymbol a$$

PCAの射影を単変量と2変量へ戻す

単変量なら、方向は $a=1$ または $-1$ しかなく、分散は $\sigma^2$ のままです。2変量なら $Z=a_1X_1+a_2X_2$ で、分散は $a_1^2\sigma_1^2+a_2^2\sigma_2^2+2a_1a_2\sigma_{12}$ です。PCAは、$a_1^2+a_2^2=1$ を保ちながら、この分散が最大になる重みを探しています。

です。$\|\boldsymbol a\|=1$ の下でこの分散を最大にする方向は、共分散行列 $\boldsymbol\Sigma$ の最大固有値に対応する固有ベクトルです。これが第一主成分です。

PCAは「見た目を2次元にする方法」ではなく、線形結合の分散を最大にする座標変換です。 特徴量の単位が大きく異なる場合、共分散行列を使うか相関行列を使うかで結果が変わります。

マハラノビス距離は分布外入力を見つける

訓練データの平均を $\boldsymbol\mu$、共分散行列を $\boldsymbol\Sigma$ とすると、

$$D_M^2(\boldsymbol x) =(\boldsymbol x-\boldsymbol\mu)^T \boldsymbol\Sigma^{-1} (\boldsymbol x-\boldsymbol\mu)$$

は、相関と尺度を補正した距離です。 化合物の分子記述子や患者のバイオマーカーが訓練集団から遠い場合、モデルが高い確率を出していても、その予測は外挿かもしれません。

これはout-of-distribution detectionとモデルの適用領域の問題です。 ただし、高次元で標本数が少ないと共分散行列が不安定になるため、正則化、次元削減、外部検証を組み合わせます。

ガウス混合モデルは所属を確率で表す

潜在クラス $Z\in\{1,\ldots,K\}$ を導入すると、

$$p(\boldsymbol x) =\sum_{k=1}^K\pi_k \phi_p(\boldsymbol x;\boldsymbol\mu_k,\boldsymbol\Sigma_k)$$

ガウス混合モデルを単変量へ戻す

単変量なら、$\phi_p$ は平均 $\mu_k$、分散 $\sigma_k^2$ の普通の正規密度です。多変量では $\mu_k$ が各特徴量の中心を並べたベクトル、$\Sigma_k$ が成分 $k$ 内の分散と相関を表します。各密度を混合比 $\pi_k$ で重み付けして足す操作自体は単変量と同じです。

です。観測後の成分所属確率は、

$$P(Z=k\mid\boldsymbol X=\boldsymbol x) =\frac{\pi_k\phi_p(\boldsymbol x;\boldsymbol\mu_k,\boldsymbol\Sigma_k)} {\sum_j\pi_j\phi_p(\boldsymbol x;\boldsymbol\mu_j,\boldsymbol\Sigma_j)}$$

となります。各標本を1つの群へ即座に固定せず、所属の曖昧さを残すソフトクラスタリングです。

混合成分が生物学的な細胞型や反応群に必ず対応するわけではありません。初期値、成分数、バッチ効果、外部マーカーによる検証が必要です。

埋め込み空間の距離は自動的に意味を持たない

深層学習では、画像、分子グラフ、遺伝子発現を低次元ベクトルへ変換した潜在表現を使います。 しかし、その空間のユークリッド距離が生物学的類似度を表すには、学習目的、尺度、訓練データが適切である必要があります。

同じ距離でも、分散が小さい方向のずれは異常性が高く、分散が大きい方向のずれは珍しくない場合があります。 マハラノビス距離やwhiteningは、この方向ごとのばらつきを補正する考え方です。

データ分割でも階層を守る

同一患者の複数時点、同一マウスの複数切片、同一ドナー由来の多数細胞をランダムに訓練・テストへ分けると、同じ上位単位の情報が両方へ入ります。 これは情報漏洩であり、汎化性能を過大評価します。

独立単位が患者なら患者単位、マウスならマウス単位、ドナーならドナー単位で分割します。 階層モデルで推定時の依存を扱うことと、group splitで評価時の漏洩を防ぐことは、別々に必要です。

16. 統計検定1級で優先する論点

1. 同時分布の台を図示し、正規化する
2. 周辺分布は不要な変数について和・積分を取る
3. 条件付き分布は同時分布を周辺分布で割る
4. 反復期待値と全分散の公式を使えるようにする
5. 共分散を $E[XY]-E[X]E[Y]$ で計算する
6. 独立と無相関を区別する
7. 和の分布を畳み込みまたは母関数で求める
8. 多次元変換では台、逆変換、ヤコビアンをセットで書く
9. 多項分布の負の共分散を指示変数から導く
10. 多変量正規分布の線形結合と条件付き分布を使う
11. 混合分布の分散を群内と群間に分解する
12. 階層を無視した疑似反復を見抜く

17. 数理統計問題

ここでは、同時分布、条件付き期待値、畳み込み、変数変換、多項分布、多変量正規分布を、抽象的な確率変数から導く練習をします。

問題1：同時確率表

同時確率が次で与えられます。

	$Y=0$	$Y=1$
$X=0$	0.42	0.08
$X=1$	0.18	0.32

次を求めてください。

1. $P(X=1)$ と $P(Y=1)$
2. $P(Y=1\mid X=1)$
3. $X,Y$ は独立か
4. $\mathrm{Cov}(X,Y)$ と相関係数

解答1

周辺確率は行・列を足して、

$$P(X=1)=0.18+0.32=0.50$$ $$P(Y=1)=0.08+0.32=0.40$$

です。条件付き確率は、

$$P(Y=1\mid X=1) =\frac{0.32}{0.50} =0.64$$

です。独立なら $P(X=1,Y=1)=0.50\times0.40=0.20$ ですが、実際は0.32なので独立ではありません。

二値変数なので、

$$E[XY]=P(X=1,Y=1)=0.32$$

です。したがって、

$$\mathrm{Cov}(X,Y) =0.32-0.50\times0.40 =0.12$$

です。また、

$$\mathrm{Var}(X)=0.5(1-0.5)=0.25$$ $$\mathrm{Var}(Y)=0.4(1-0.4)=0.24$$

なので、

$$\mathrm{Corr}(X,Y) =\frac{0.12}{\sqrt{0.25\times0.24}} \approx0.490$$

です。

問題2：三角形の台

$$f_{X,Y}(x,y)=c \qquad(0<x<y<1)$$

が同時密度となるように $c$ を定め、周辺密度 $f_X(x),f_Y(y)$ と条件付き密度 $f_{X\mid Y}(x\mid y)$ を求めてください。

解答2

正規化より、

$$1=\int_0^1\int_0^y c\,dx\,dy =\frac c2$$

なので $c=2$ です。

$$f_X(x)=\int_x^1 2\,dy=2(1-x) \qquad(0<x<1)$$ $$f_Y(y)=\int_0^y2\,dx=2y \qquad(0<y<1)$$

です。したがって、

$$f_{X\mid Y}(x\mid y) =\frac{2}{2y} =\frac1y \qquad(0<x<y)$$

です。

問題3：条件付き期待値

問題1の表について、$E[Y\mid X]$ と $E\{E[Y\mid X]\}$ を求めてください。

解答3

$$E[Y\mid X=0] =P(Y=1\mid X=0) =\frac{0.08}{0.50} =0.16$$ $$E[Y\mid X=1]=0.64$$

です。したがって、

$$E\{E[Y\mid X]\} =0.50\times0.16+0.50\times0.64 =0.40$$

となり、$E[Y]=P(Y=1)=0.40$ と一致します。

問題4：全分散の公式

問題1の表について、全分散の公式を使って $\mathrm{Var}(Y)$ を求めてください。

解答4

条件付き分散は、

$$\mathrm{Var}(Y\mid X=0) =0.16(1-0.16) =0.1344$$ $$\mathrm{Var}(Y\mid X=1) =0.64(1-0.64) =0.2304$$

です。よって、

$$E\{\mathrm{Var}(Y\mid X)\} =0.5(0.1344)+0.5(0.2304) =0.1824$$

です。条件付き平均は0.16と0.64で、全体平均は0.40なので、

$$\mathrm{Var}\{E[Y\mid X]\} =0.5(0.16-0.40)^2 +0.5(0.64-0.40)^2 =0.0576$$

です。したがって、

$$\mathrm{Var}(Y)=0.1824+0.0576=0.24$$

となり、ベルヌーイ分散 $0.4(1-0.4)$ と一致します。

問題5：線形結合の平均と分散

$$E[X]=2,\quad E[Y]=3,$$ $$\mathrm{Var}(X)=4,\quad \mathrm{Var}(Y)=9,\quad \mathrm{Cov}(X,Y)=3$$

とします。$S=2X-Y$ の平均と分散を求めてください。

解答5

$$E[S]=2E[X]-E[Y]=4-3=1$$

です。分散は、

$$\begin{aligned} \mathrm{Var}(2X-Y) &=4\mathrm{Var}(X) +\mathrm{Var}(Y) -4\mathrm{Cov}(X,Y)\\ &=4\times4+9-4\times3\\ &=13 \end{aligned}$$

です。係数の積 $2\times(-1)$ を2倍して、共分散項が $-4\mathrm{Cov}(X,Y)$ となる点に注意します。

問題6：離散畳み込み

独立な確率変数 $X,Y$ が、

$$P(X=0)=P(Y=0)=\frac12, \qquad P(X=1)=P(Y=1)=\frac12$$

を満たします。$S=X+Y$ の分布を求めてください。

解答6

$$P(S=0)=P(X=0,Y=0)=\frac14$$ $$P(S=1) =P(0,1)+P(1,0) =\frac14+\frac14 =\frac12$$ $$P(S=2)=P(1,1)=\frac14$$

です。したがって $S\sim\mathrm{Bin}(2,1/2)$ です。

問題7：一様分布の和

$X,Y$ は独立に $U(0,1)$ に従います。$S=X+Y$ の密度を求めてください。

解答7

畳み込みより、

$$f_S(s)=\int f_X(x)f_Y(s-x)\,dx$$

です。積分される値が1になるには、

$$0<x<1, \qquad 0<s-x<1$$

が必要です。

$0<s<1$ では $0<x<s$ なので、

$$f_S(s)=\int_0^s1\,dx=s$$

です。$1\leq s<2$ では $s-1<x<1$ なので、

$$f_S(s)=\int_{s-1}^1 1\,dx=2-s$$

です。したがって、

$$f_S(s) = \begin{cases} s & (0<s<1)\\ 2-s & (1\leq s<2)\\ 0 & (\text{otherwise}) \end{cases}$$

です。

問題8：指数分布の和

$X,Y$ は独立に率 $\lambda$ の指数分布に従います。$S=X+Y$ の密度を畳み込みで求めてください。

解答8

$s>0$ で、

$$\begin{aligned} f_S(s) &=\int_0^s \lambda e^{-\lambda x} \lambda e^{-\lambda(s-x)}\,dx\\ &=\lambda^2e^{-\lambda s} \int_0^s1\,dx\\ &=\lambda^2s e^{-\lambda s} \end{aligned}$$

です。これは形状2、率 $\lambda$ のガンマ分布です。

問題9：和と比への変換

$X,Y$ は独立に率1の指数分布に従います。

$$S=X+Y,\qquad U=\frac{X}{X+Y}$$

とします。$(S,U)$ の同時密度を求め、$S,U$ が独立であることを示してください。

解答9

逆変換は、

$$X=SU,\qquad Y=S(1-U)$$

です。台は、

$$s>0,\qquad0<u<1$$

です。ヤコビアンは、

$$\left| \frac{\partial(x,y)}{\partial(s,u)} \right| = \left| \begin{matrix} u&s\\ 1-u&-s \end{matrix} \right| =s$$

です。したがって、

$$\begin{aligned} f_{S,U}(s,u) &=e^{-su}e^{-s(1-u)}s\\ &=se^{-s} \end{aligned}$$

です。これは、

$$\{se^{-s}\}\times1$$

と分解できるので、

$$S\sim\mathrm{Gamma}(2,\text{rate }1), \qquad U\sim U(0,1)$$

であり、$S,U$ は独立です。

問題10：Box–Muller変換

$$U_1=e^{-1}, \qquad U_2=\frac18$$

について、次を求めてください。

1. 半径 $R$ と角度 $\Theta$
2. Box–Muller変換後の $Z_1,Z_2$
3. $Z_1^2+Z_2^2=-2\log U_1$ が成り立つことの確認
4. $U_1$ と $U_2$ がそれぞれ変換後の点の何を決めているか
5. 極座標から直交座標へのヤコビアンの絶対値
6. $(U_1,U_2)$ から $(Z_1,Z_2)$ への変換全体のヤコビアンと、その逆数

解答10

半径は、

$$R=\sqrt{-2\log(e^{-1})}=\sqrt2$$

です。角度は、

$$\Theta=2\pi\times\frac18=\frac{\pi}{4}$$

です。したがって、

$$Z_1=\sqrt2\cos\frac{\pi}{4}=1$$ $$Z_2=\sqrt2\sin\frac{\pi}{4}=1$$

です。

平方和を計算すると、

$$Z_1^2+Z_2^2=1^2+1^2=2$$

です。一方、

$$-2\log U_1 =-2\log(e^{-1}) =2$$

なので、

$$Z_1^2+Z_2^2=-2\log U_1$$

を確認できます。

$U_1$ は、

$$R=\sqrt{-2\log U_1}$$

を通して原点からの距離を決めます。$U_1$ が0へ近づくほど $R$ は大きくなります。 $U_2$ は、

$$\Theta=2\pi U_2$$

を通して方向を決め、半径は変えません。

極座標から直交座標へのヤコビアンは、

$$\left| \frac{\partial(z_1,z_2)} {\partial(r,\theta)} \right| =r$$

です。この問題では $R=\sqrt2$ なので、

$$\left| \frac{\partial(z_1,z_2)} {\partial(r,\theta)} \right| =\sqrt2$$

です。これは $(R,\Theta)$ 空間の微小長方形 $dR,d\Theta$ が、$(Z_1,Z_2)$ 平面で面積約 $\sqrt2,dR,d\Theta$ の扇形へ変わることを表します。

変換全体のヤコビアンは、

$$\left| \frac{\partial(z_1,z_2)} {\partial(u_1,u_2)} \right| =\frac{2\pi}{u_1}$$

でした。$u_1=e^{-1}$ を代入すると、

$$\left| \frac{\partial(z_1,z_2)} {\partial(u_1,u_2)} \right| =2\pi e$$

です。密度変換で $(Z_1,Z_2)$ の密度を求めるときに使う逆向きのヤコビアンは、

$$\left| \frac{\partial(u_1,u_2)} {\partial(z_1,z_2)} \right| =\frac{1}{2\pi e}$$

です。

問題11：多項分布の確率

細胞の表現型がA、B、Cに分類される確率を、

$$(p_A,p_B,p_C)=(0.5,0.3,0.2)$$

とします。独立に10細胞を観察したとき、個数が $(5,3,2)$ となる確率を求めてください。

解答11

$$P(N_A=5,N_B=3,N_C=2) = \frac{10!}{5!3!2!} (0.5)^5(0.3)^3(0.2)^2$$

です。多項係数は、

$$\frac{10!}{5!3!2!}=2520$$

なので、

$$2520(0.5)^5(0.3)^3(0.2)^2 \approx0.0851$$

です。

問題12：多項分布の共分散

問題11で、$E[N_A]$、$\mathrm{Var}(N_A)$、$\mathrm{Cov}(N_A,N_B)$ を求めてください。

解答12

$$E[N_A]=np_A=10\times0.5=5$$ $$\mathrm{Var}(N_A) =np_A(1-p_A) =10\times0.5\times0.5 =2.5$$ $$\mathrm{Cov}(N_A,N_B) =-np_Ap_B =-10\times0.5\times0.3 =-1.5$$

です。合計細胞数が固定されているため、カテゴリ間の共分散は負です。

問題13：多変量正規分布の線形結合

$$\begin{pmatrix}X\\Y\end{pmatrix} \sim N_2\left[ \begin{pmatrix}10\\20\end{pmatrix}, \begin{pmatrix}4&3\\3&9\end{pmatrix} \right]$$

とします。$S=X+2Y$ の分布を求めてください。

解答13

平均は、

$$E[S]=10+2\times20=50$$

です。$\mathbf a=(1,2)^\mathsf T$ とすると、

$$\begin{aligned} \mathrm{Var}(S) &=\mathbf a^\mathsf T \boldsymbol\Sigma\mathbf a\\ &=1^2\times4 +2^2\times9 +2\times1\times2\times3\\ &=52 \end{aligned}$$

です。したがって、

$$S\sim N(50,52)$$

です。

問題14：二変量正規の条件付き分布

$X,Y$ は平均0、分散1、相関係数 $\rho=0.8$ の二変量正規分布に従います。 $X=1.5$ のときの $Y$ の条件付き分布を求めてください。

解答14

$$E[Y\mid X=1.5] =\rho x =0.8\times1.5 =1.2$$ $$\mathrm{Var}(Y\mid X=1.5) =1-\rho^2 =1-0.64 =0.36$$

なので、

$$Y\mid X=1.5\sim N(1.2,0.36)$$

です。

問題15：2成分混合分布

確率0.7で $N(0,1)$、確率0.3で $N(4,1)$ から値を生成します。 混合分布の平均と分散を求めてください。

解答15

平均は、

$$\mu=0.7\times0+0.3\times4=1.2$$

です。群内分散の平均は、

$$0.7\times1+0.3\times1=1$$

です。群間平均による分散は、

$$0.7(0-1.2)^2+0.3(4-1.2)^2 =3.36$$

です。したがって、

$$\mathrm{Var}(X)=1+3.36=4.36$$

です。

問題16：同じユークリッド距離、異なるマハラノビス距離

平均ベクトルが $\boldsymbol{0}$、各変数の分散が1、相関係数が $\rho=0.8$ の二変量正規分布を考えます。 次の2点について、原点からのユークリッド距離とマハラノビス距離の2乗を求めてください。

$$\boldsymbol{x}_A=(2,2)^\top, \qquad \boldsymbol{x}_B=(2,-2)^\top$$

さらに、自由度2のカイ二乗分布の上側5%点 $5.991$ と比較し、どちらが95%確率楕円の外側にあるか説明してください。

解答16

どちらの点もユークリッド距離は、

$$\|\boldsymbol{x}_A\|_2 =\|\boldsymbol{x}_B\|_2 =\sqrt{2^2+2^2} =\sqrt{8}$$

であり、原点からの直線距離だけでは区別できません。

各変数の標準偏差が1なので、マハラノビス距離の2乗は、

$$D^2 =\frac{x_1^2-2\rho x_1x_2+x_2^2}{1-\rho^2}$$

です。$\boldsymbol{x}_A=(2,2)^\top$ では、

$$\begin{aligned} D_A^2 &=\frac{2^2-2(0.8)(2)(2)+2^2}{1-0.8^2}\\ &=\frac{8-6.4}{0.36}\\ &=\frac{1.6}{0.36}\\ &\approx4.44 \end{aligned}$$

です。一方、$\boldsymbol{x}_B=(2,-2)^\top$ では、

$$\begin{aligned} D_B^2 &=\frac{2^2-2(0.8)(2)(-2)+(-2)^2}{1-0.8^2}\\ &=\frac{8+6.4}{0.36}\\ &=\frac{14.4}{0.36}\\ &=40 \end{aligned}$$

です。したがって、

$$D_A^2=4.44<5.991, \qquad D_B^2=40>5.991$$

より、$\boldsymbol{x}_A$ は95%確率楕円の内側、$\boldsymbol{x}_B$ は外側です。 正の相関が強い集団では「両方とも高い」は起こりやすい一方、「片方が高く、もう片方が低い」は非常に起こりにくいためです。

問題17：階層モデルの分散、ICC、部分プーリング

次の正規階層モデルを考えます。

$$Y_{ij}=\mu+u_i+\varepsilon_{ij}$$ $$u_i\sim N(0,0.64), \qquad \varepsilon_{ij}\sim N(0,1.44)$$

各群の標本数は $n_i=4$、全体平均は $\mu=1$ とします。ある群の標本平均が $\bar Y_i=3$ でした。

1. $Y_{ij}$ の分散
2. 同じ群内の2観測の共分散
3. ICC
4. 部分プーリング後の群平均

を求めてください。

解答17

群間分散は $\tau^2=0.64$、群内分散は $\sigma^2=1.44$ です。

観測値の分散は、

$$\begin{aligned} \mathrm{Var}(Y_{ij}) &=\mathrm{Var}(u_i)+\mathrm{Var}(\varepsilon_{ij})\\ &=0.64+1.44\\ &=2.08 \end{aligned}$$

です。

同じ群の2観測は $u_i$ を共有するため、

$$\mathrm{Cov}(Y_{ij},Y_{ik}) =\mathrm{Var}(u_i) =0.64$$

です。したがってICCは、

$$\begin{aligned} \mathrm{ICC} &=\frac{\tau^2}{\tau^2+\sigma^2}\\ &=\frac{0.64}{0.64+1.44}\\ &=\frac{0.64}{2.08}\\ &\approx0.308 \end{aligned}$$

です。

部分プーリングで群自身の平均に掛かる重みは、

$$\begin{aligned} B_i &=\frac{n_i\tau^2}{n_i\tau^2+\sigma^2}\\ &=\frac{4(0.64)}{4(0.64)+1.44}\\ &=\frac{2.56}{4.00}\\ &=0.64 \end{aligned}$$

です。したがって、

$$\begin{aligned} E[\theta_i\mid\bar Y_i] &=B_i\bar Y_i+(1-B_i)\mu\\ &=0.64(3)+0.36(1)\\ &=1.92+0.36\\ &=2.28 \end{aligned}$$

です。生の群平均3をそのまま採用せず、標本数と分散成分に応じて全体平均1の方向へ縮めています。

問題18：Cauchy—Schwarzの不等式を二次式から導く

$E[X^2]<\infty$、$E[Y^2]<\infty$ とします。二次式

$$h(t)=E[\{(X-E[X])-t(Y-E[Y])\}^2]$$

を用いて、

$$|\operatorname{Cov}(X,Y)| \leq\sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}$$

を示し、等号成立条件を答えてください。

解答18

$U=X-E[X]$、$V=Y-E[Y]$ とおくと、

$$\begin{aligned} h(t) &=E[(U-tV)^2]\\ &=\operatorname{Var}(X)-2t\operatorname{Cov}(X,Y) +t^2\operatorname{Var}(Y). \end{aligned}$$

$h(t)$ は平方の期待値なので、すべての $t$ に対して非負です。$\operatorname{Var}(Y)>0$ のとき、平方完成すると、

$$\begin{aligned} h(t) &=\operatorname{Var}(Y) \left(t-\frac{\operatorname{Cov}(X,Y)}{\operatorname{Var}(Y)}\right)^2\\ &\quad+\operatorname{Var}(X) -\frac{\operatorname{Cov}(X,Y)^2}{\operatorname{Var}(Y)}. \end{aligned}$$

最小値も非負なので、

$$\operatorname{Cov}(X,Y)^2 \leq\operatorname{Var}(X)\operatorname{Var}(Y).$$

両辺の平方根を取れば目的の不等式を得ます。

等号が成立するのは最小値が0、すなわち、ある定数 $c$ が存在して、

$$X-E[X]=c\{Y-E[Y]\}\quad\text{a.s.}$$

となるときです。どちらかの分散が0の場合も等号になります。

問題19：条件付き期待値は残差と直交する

$X,Y$ は二乗可積分とし、

$$R=Y-E[Y\mid X]$$

とおきます。次を示してください。

1. $\operatorname{Cov}(X,R)=0$
2. $\operatorname{Var}(R)=E[\operatorname{Var}(Y\mid X)]$

解答19

まず、条件付き期待値の定義から、

$$E[R\mid X] =E[Y\mid X]-E[Y\mid X] =0.$$

したがって、反復期待値の法則により、

$$E[R]=E[E[R\mid X]]=0$$

です。また、$X$ は $X$ の関数なので条件の外へ出せます。

$$\begin{aligned} E[XR] &=E[E[XR\mid X]]\\ &=E[XE[R\mid X]]\\ &=0. \end{aligned}$$

よって、

$$\operatorname{Cov}(X,R) =E[XR]-E[X]E[R] =0.$$

次に、$E[R\mid X]=0$ なので、

$$\operatorname{Var}(R) =E[R^2].$$

条件付きで見ると $E[Y\mid X]$ は定数であるため、

$$E[R^2\mid X] =E[\{Y-E[Y\mid X]\}^2\mid X] =\operatorname{Var}(Y\mid X).$$

再び期待値を取れば、

$$\operatorname{Var}(R) =E[\operatorname{Var}(Y\mid X)]$$

です。これは回帰における「説明できなかったばらつき」に対応します。

問題20：カイ二乗分布の和と比を同時変換する

独立な確率変数

$$A\sim\chi_m^2, \qquad B\sim\chi_n^2$$

に対して、

$$S=A+B, \qquad R=\frac{A}{A+B}$$

とおきます。$(S,R)$ の同時密度を求め、$S$ と $R$ の分布および独立性を示してください。

解答20

逆変換は、

$$A=RS, \qquad B=(1-R)S$$

であり、台は $s>0$、$0<r<1$ です。Jacobianは、

$$\left| \frac{\partial(a,b)}{\partial(s,r)} \right| =\left| \begin{matrix} r&s\\ 1-r&-s \end{matrix} \right| =s.$$

独立性とカイ二乗密度を使うと、

$$\begin{aligned} f_{S,R}(s,r) &=\frac{(rs)^{m/2-1}\{(1-r)s\}^{n/2-1}e^{-s/2}s} {2^{(m+n)/2}\Gamma(m/2)\Gamma(n/2)}\\ &=\left\{ \frac{s^{(m+n)/2-1}e^{-s/2}} {2^{(m+n)/2}\Gamma((m+n)/2)} \right\}\\ &\quad\times\left\{ \frac{\Gamma((m+n)/2)}{\Gamma(m/2)\Gamma(n/2)} r^{m/2-1}(1-r)^{n/2-1} \right\}. \end{aligned}$$

同時密度が $s$ だけの密度と $r$ だけの密度の積に分解できたので、

$$S\sim\chi_{m+n}^2, \qquad R\sim\operatorname{Beta}\left(\frac m2,\frac n2\right), \qquad S\perp R$$

です。

問題21：全共分散の公式と潜在変数による相関

1. 次の全共分散の公式を示してください。

$$\operatorname{Cov}(X,Y) =E[\operatorname{Cov}(X,Y\mid Z)] +\operatorname{Cov}(E[X\mid Z],E[Y\mid Z])$$

2. $Z\sim N(0,\tau^2)$ とし、$Z$ を与えたとき $X,Y$ は独立に、

$$X\mid Z\sim N(Z,\sigma^2), \qquad Y\mid Z\sim N(Z,\sigma^2)$$

に従うとします。$\operatorname{Cov}(X,Y)$ と相関係数を求めてください。

解答21

条件付き共分散の定義から、

$$E[XY\mid Z] =\operatorname{Cov}(X,Y\mid Z) +E[X\mid Z]E[Y\mid Z]$$

です。両辺の期待値を取り、$E[XY]-E[X]E[Y]$ を作ると、

$$\begin{aligned} \operatorname{Cov}(X,Y) &=E[\operatorname{Cov}(X,Y\mid Z)]\\ &\quad+E[E[X\mid Z]E[Y\mid Z]]-E[X]E[Y]\\ &=E[\operatorname{Cov}(X,Y\mid Z)]\\ &\quad+\operatorname{Cov}(E[X\mid Z],E[Y\mid Z]). \end{aligned}$$

階層モデルでは、条件付き独立性より、

$$\operatorname{Cov}(X,Y\mid Z)=0$$

です。一方、$E[X\mid Z]=E[Y\mid Z]=Z$ なので、

$$\operatorname{Cov}(X,Y)=\operatorname{Var}(Z)=\tau^2.$$

また、全分散の公式より、

$$\operatorname{Var}(X) =E[\sigma^2]+\operatorname{Var}(Z) =\sigma^2+\tau^2$$

で、$Y$ も同じです。したがって、

$$\operatorname{Corr}(X,Y) =\frac{\tau^2}{\sigma^2+\tau^2}.$$

条件付きでは独立でも、共通の潜在変数を周辺化すると正の相関が生じます。

発展問題22：Dirichlet分布の周辺・条件付き分布

$(P_1,P_2,P_3)\sim\operatorname{Dirichlet}(a,b,c)$ とし、$P_1+P_2+P_3=1$ とします。

1. $P_1$ の周辺分布を求めてください。
2. $U=P_2/(1-P_1)$ の分布を求め、$P_1$ と $U$ が独立であることを示してください。
3. $\operatorname{Cov}(P_1,P_2)$ を求めてください。

解答22

$x=P_1$、$u=P_2/(1-P_1)$ とおくと、

$$P_2=(1-x)u, \qquad P_3=(1-x)(1-u)$$

です。台は $0<x<1$、$0<u<1$ で、$(P_1,P_2)$ から $(x,u)$ への逆変換のJacobianは、

$$\left| \frac{\partial(P_1,P_2)}{\partial(x,u)} \right|=1-x$$

です。Dirichlet密度に代入すると、

$$\begin{aligned} f_{P_1,U}(x,u) &=\frac{\Gamma(a+b+c)}{\Gamma(a)\Gamma(b)\Gamma(c)} x^{a-1}\{(1-x)u\}^{b-1}\\ &\quad\times\{(1-x)(1-u)\}^{c-1}(1-x)\\ &=\left\{ \frac{x^{a-1}(1-x)^{b+c-1}}{B(a,b+c)} \right\} \left\{ \frac{u^{b-1}(1-u)^{c-1}}{B(b,c)} \right\}. \end{aligned}$$

したがって、

$$P_1\sim\operatorname{Beta}(a,b+c), \qquad U\sim\operatorname{Beta}(b,c), \qquad P_1\perp U.$$

$\alpha_0=a+b+c$ とおくと、Dirichlet分布の積率は、

$$E[P_1]=\frac{a}{\alpha_0}, \qquad E[P_2]=\frac{b}{\alpha_0}, \qquad E[P_1P_2]=\frac{ab}{\alpha_0(\alpha_0+1)}$$

なので、

$$\operatorname{Cov}(P_1,P_2) =-\frac{ab}{\alpha_0^2(\alpha_0+1)}.$$

成分和が1に固定されるため、1成分が増えると他成分に負の共分散が生じます。

18. 応用統計問題（医薬・生命科学）

ここでは、バイオマーカー、有害事象、細胞・個体の階層構造など、医薬・生命科学データの実験単位と不確実性をモデルから読み取ります。

問題1：バイオマーカー層別と全分散

患者の30%がバイオマーカー陽性です。 陽性群の薬効スコアの平均・分散はそれぞれ8、4、陰性群ではそれぞれ5、9とします。 全患者での平均と分散を求めてください。

解答1

陽性を $Z=1$ とすると、

$$E[Y] =0.3\times8+0.7\times5 =5.9$$

です。群内分散の平均は、

$$E\{\mathrm{Var}(Y\mid Z)\} =0.3\times4+0.7\times9 =7.5$$

です。群間平均の分散は、

$$\begin{aligned} \mathrm{Var}\{E[Y\mid Z]\} &=0.3(8-5.9)^2+0.7(5-5.9)^2\\ &=1.89 \end{aligned}$$

です。したがって、

$$\mathrm{Var}(Y)=7.5+1.89=9.39$$

です。

問題2：有害事象グレードの多項分布

100人の患者について、有害事象なし、軽度、中等度、重度の確率が、

$$(0.55,0.25,0.15,0.05)$$

とします。

1. 重度の人数 $N_4$ の平均と分散
2. 中等度以上の人数 $N_3+N_4$ の分布
3. 軽度と重度の人数の共分散

を求めてください。

解答2

重度の人数の周辺分布は、

$$N_4\sim\mathrm{Bin}(100,0.05)$$

なので、

$$E[N_4]=5, \qquad \mathrm{Var}(N_4)=100(0.05)(0.95)=4.75$$

です。

中等度以上の確率は、

$$0.15+0.05=0.20$$

なので、

$$N_3+N_4\sim\mathrm{Bin}(100,0.20)$$

です。軽度と重度の共分散は、

$$\mathrm{Cov}(N_2,N_4) =-100(0.25)(0.05) =-1.25$$

です。

問題3：階層構造と疑似反復

薬剤群3匹、対照群3匹のマウスから、それぞれ100細胞を測定しました。 600細胞を独立標本として2群比較する解析の問題点を説明し、適切な解析単位とモデル案を示してください。

解答3

同じマウス由来の100細胞は、遺伝背景、処置、組織環境、標本作製条件を共有するため独立ではありません。 600細胞を独立標本とみなすと、実際の独立単位より標本数を過大に数え、標準誤差を過小評価する疑似反復になります。

処置がマウスに割り付けられているなら、処置効果の独立な反復単位は基本的にマウスです。 解析案としては、

• マウスごとに細胞測定値を要約し、$n=3$ 対 $n=3$ として比較する
• 細胞をレベル1、マウスをレベル2とする混合効果モデルを用いる
• 細胞型やバッチがあるなら、研究目的とデザインに応じて追加階層を組み込む

ことが考えられます。 細胞数を増やすことは各マウス内の推定精度を高めますが、マウス数を増やすことと同じではありません。

問題4：ベータ二項分布の過分散

施設ごとの有害事象発現確率が、

$$P\sim\mathrm{Beta}(4,16)$$

に従い、各施設で $n=20$ 人を観察するとします。 施設ごとの発現人数 $Y$ の平均と分散を求め、成功確率を固定した二項分布と比較してください。

解答4

$$\mu=\frac{4}{4+16}=0.2, \qquad \kappa=20$$

です。平均は、

$$E[Y]=n\mu=20\times0.2=4$$

です。ベータ二項分布の分散は、

$$\begin{aligned} \mathrm{Var}(Y) &=n\mu(1-\mu) \times\frac{\kappa+n}{\kappa+1}\\ &=20(0.2)(0.8) \times\frac{20+20}{20+1}\\ &\approx6.095 \end{aligned}$$

です。一方、成功確率を0.2に固定した二項分布の分散は、

$$20(0.2)(0.8)=3.2$$

です。施設ごとの発現確率の違いにより、分散が約1.90倍になっています。

問題5：多変量バイオマーカーの条件付き予測

標準化した2つのバイオマーカー $X,Y$ が相関0.6の二変量正規分布に従うとします。

1. $X=2$ の患者における $Y$ の条件付き平均と分散
2. 条件付け前後で分散が何%減少するか

を求めてください。

解答5

$$Y\mid X=2 \sim N(0.6\times2,1-0.6^2) =N(1.2,0.64)$$

です。したがって条件付き平均は1.2、分散は0.64です。 条件付け前の分散は1なので、分散の減少率は、

$$\frac{1-0.64}{1}\times100=36\%$$

です。 これは相関があるとき、$X$ の情報が $Y$ の予測不確実性を減らすことを表します。 ただし、予測精度の評価には、正規性、線形性、外部検証、測定誤差も確認する必要があります。

問題6：AUCとCmaxの組み合わせ異常

ある薬物動態試験で、用量補正後の $\log(\mathrm{AUC})$ と $\log(C_{\max})$ を標準化した値をそれぞれ $X,Y$ とします。 基準集団では平均がともに0、分散がともに1、相関係数が $\rho=0.7$ でした。 ある被験者の観測値が、

$$(X,Y)=(1.5,-1.0)$$

であったとします。

1. この被験者のマハラノビス距離の2乗 $D^2$ を求めてください。
2. 二変量正規分布を仮定し、$\chi_2^2$ の上側5%点 $5.991$ と比較してください。
3. 実データで確認すべき項目を3つ以上挙げてください。

解答6

標準化済みなので、

$$D^2 =\frac{X^2-2\rho XY+Y^2}{1-\rho^2}$$

を使います。値を代入すると、

$$\begin{aligned} D^2 &=\frac{1.5^2-2(0.7)(1.5)(-1.0)+(-1.0)^2}{1-0.7^2}\\ &=\frac{2.25+2.10+1.00}{0.51}\\ &=\frac{5.35}{0.51}\\ &\approx10.49 \end{aligned}$$

したがって、

$$D\approx\sqrt{10.49}\approx3.24$$

です。また、

$$D^2=10.49>5.991$$

なので、この点は基準集団の95%確率楕円の外側にあります。 AUCが高い一方で $C_{\max}$ が低いという組み合わせは、正の相関をもつ基準集団では起こりにくいと解釈できます。

ただし、これだけで「異常な薬物動態」や「疾患による変化」と断定してはいけません。少なくとも次を確認します。

• 投与量、服薬遵守、投与・採血時刻の記録
• 単位、転記、欠測補完、用量補正などのデータ処理
• 測定法の定量下限、精度管理、バッチ効果
• 吸収速度、製剤、食事、併用薬などの臨床条件
• 腎機能、肝機能、体重、遺伝型などの患者背景
• 基準集団の選び方、正規性、共分散行列の推定精度

マハラノビス距離は、確認すべき症例を見つけるスクリーニング指標です。原因を特定する指標ではないため、記録・測定・生物学的背景を順に検証します。

問題7：混合分布か階層モデルか

次の研究場面で、混合分布、階層モデル、または両方のどれを考えるべきか答え、理由を説明してください。

1. 由来不明の細胞1000個について薬剤応答量だけがあり、感受性群と耐性群の混在を調べたい
2. 6匹のマウスから各500細胞を測定し、マウスIDが記録されている
3. 10施設の患者を追跡し、施設内にもレスポンダーと非レスポンダーがいる可能性がある

解答7

1では、各細胞の表現型ラベルが観測されておらず、未知の感受性群・耐性群を推定したいので混合分布が候補です。ただし、二峰性だけで生物学的表現型を断定せず、マーカーや機能試験で検証します。

2では、細胞がどのマウスに由来するか既知であり、細胞がマウス内に入れ子になっているので階層モデルが候補です。500細胞を500匹の独立マウスとして扱ってはいけません。処置効果の一般化対象がマウスなら、独立単位は基本的に6匹です。

3では、患者が施設内に入れ子になっているため施設差を扱う階層モデルが必要です。さらに、施設内で所属不明のレスポンダー・非レスポンダーが混在すると仮定するなら、潜在クラスをもつ混合分布も組み合わせることになります。

実際のモデル選択では、次を順に確認します。

• 所属IDが観測されている階層は何か
• 未観測のサブタイプを仮定する根拠があるか
• 独立な実験単位は患者、個体、施設のどれか
• 技術反復と生物学的反復を区別できているか
• 複雑なモデルを支える標本数と検証データがあるか

問題8：毒性予測と適用領域

標準化した2つの化学記述子 $X_1,X_2$ の学習データで、平均がともに0、分散がともに1、相関係数が $\rho=0.65$ だったとします。

候補化合物の記述子が、

$$(X_1,X_2)=(1.5,-1.0)$$

であり、説明用の仮想毒性モデルが、

$$\operatorname{logit} \{P(\text{毒性})\} =-0.3+0.9X_1-0.6X_2$$

であったとします。

1. マハラノビス距離の2乗 $D^2$ を求めてください。
2. $\chi_2^2$ の95%点 $5.991$ と比較してください。
3. 仮想モデルの毒性確率を求めてください。
4. この候補化合物について、どのように報告し、何を追加確認すべきか説明してください。

解答8

マハラノビス距離の2乗は、

$$D^2 = \frac{x_1^2-2\rho x_1x_2+x_2^2} {1-\rho^2}$$

です。値を代入すると、

$$\begin{aligned} D^2 &= \frac{ 1.5^2-2(0.65)(1.5)(-1.0)+(-1.0)^2 }{ 1-0.65^2 }\\ &= \frac{2.25+1.95+1.00}{0.5775}\\ &= \frac{5.20}{0.5775}\\ &\approx9.00 \end{aligned}$$

です。

$$D^2\approx9.00>5.991$$

なので、二変量正規分布と母数既知を仮定した95%化学空間の外側です。

仮想毒性モデルの線形予測子は、

$$\begin{aligned} \eta &=-0.3+0.9(1.5)-0.6(-1.0)\\ &=-0.3+1.35+0.60\\ &=1.65 \end{aligned}$$

です。したがって、

$$\begin{aligned} P(\text{毒性}) &=\frac{1}{1+e^{-1.65}}\\ &\approx0.839 \end{aligned}$$

です。

報告では、

仮想モデルは毒性確率約83.9%を出力したが、候補化合物はマハラノビス距離による95%適用領域外であり、外挿予測である

と、予測値と適用領域を分けて記載します。

少なくとも次を確認します。

• 学習データ中の近傍化合物と構造類似性
• 毒性エンドポイントの定義、ラベル品質、クラス不均衡
• 構造アラート、反応性官能基、代謝活性化の可能性
• in vitro assayの濃度反応、細胞毒性、測定干渉
• 別モデル、read-across、外部データによる予測の再現性
• 実験確認の優先順位と、ヒト曝露量との関係

高い予測値を無視してはいけませんが、「83.9%を高精度に信頼できる」とも結論できません。適用領域外という情報を、追加評価を優先する警告として使います。

まとめ

この章では、複数の確率変数を同時に扱うための考え方を整理しました。

• 同時分布は複数変数の組み合わせ全体を記述する
• 周辺分布は不要な変数を足す・積分することで得る
• 条件付き分布は、得られた情報によって分布を更新する
• 反復期待値と全分散の公式は、層別・階層構造の理解につながる
• 共分散と相関は線形な連動を要約するが、独立性や因果性とは異なる
• 独立な和の分布は畳み込みで求め、母関数の積でも理解できる
• 多次元変換では逆変換、台、ヤコビアンを必ず確認する
• 多項分布はカテゴリ別個数を扱い、合計制約からカテゴリ間に負の共分散が生じる
• 多変量正規分布は線形変換、周辺化、条件付けで閉じている
• 混合分布は潜在集団、階層モデルは個体・施設・バッチなどの構造を表す
• 医薬生物データでは、観測数と独立な実験単位を区別する

多次元分布の本質は、変数を増やすことではありません。 「何を同時に観測し、何を条件として固定し、どのばらつきを分離したいか」を数式で明確にすることです。