統計学の基礎事項まとめその他の話題（理工学分野）

統計学の基礎をすぐに思い出せるような、まとめチートシートを作成しました。

「確率分布編」「推定・検定編」「線形回帰編」「その他の話題（理工学分野）」の４部構成です。本ページは「その他の話題（理工学分野）」です。

注意

あくまで個人用のメモであり、誤植や不正確な情報が含まれている可能性があります。適宜更新していきます。

A その他の話題（理工学分野）

主成分分析

多変量データでは、変数の数が多いほど全体像をつかみにくくなる。主成分分析は、元の変数を互いに直交する新しい軸に取り直し、データのばらつきが大きい方向から順に見ることで、できるだけ情報を失わずに低次元でデータを眺めるための方法である。

まずは、$D$ 次元ベクトルのデータ $(\boldsymbol x_n)_{n=1,2,\cdots,N}$ を新たな基底 $\left< \boldsymbol u _i \right>=\left< \boldsymbol u_1, \cdots, \boldsymbol u_D \right>$ で表現することを考えよう。ここで基底 $\left< \boldsymbol u_i \right>$ は正規直交基底とする。つまり

\begin{align} \boldsymbol u_i^\top \boldsymbol u_j = \delta_{ij}, \quad \sum_i \boldsymbol u_i \boldsymbol u_i^\top = 1 \end{align}

である。ただしこれまでと同様、単位行列 $I$ は誤解のない限りたんに $1$ と書く。後半の完全性の式を用いると、データ $ \boldsymbol x_n \in \mathbb{R}^D $ は基底 $\left< \boldsymbol u _i \right>$ を用いて一意に展開できる。

\begin{align} \boldsymbol x_n = \sum_i \boldsymbol u_i (\boldsymbol u_i^\top \boldsymbol x_n) \end{align}

さていま、データ $\boldsymbol x_n$ を $M$ ($\leq D$) 個の基底ベクトル $\boldsymbol u_1, \ldots, \boldsymbol u_M$ を用いて、なるべく近似したいとしよう。つまり、

\begin{align} \boldsymbol x_n \simeq \sum_{i=1}^M \boldsymbol u_i (\boldsymbol u_i^\top \boldsymbol x_n) + C \end{align}

となるような基底 $\left< \boldsymbol u_1, \cdots, \boldsymbol u_M, \cdots, \boldsymbol u_D \right>$ を見つけたいとしよう。ただし $C$ は $n$ によらない定数である。言い換えると、$C$ が定数というのは「 $i \geq M+1$ のとき、どのデータ $\boldsymbol x_n$ でも基底 $\boldsymbol u_i$ の成分がほとんど一定」ということを示している。このような場合、データ $\boldsymbol x_n$ の $\boldsymbol u_i$ 成分の分散は、ほとんどゼロに近いであろう。

このような基底はどのようにすれば作れるのだろうか？これについては、誤差関数とLagrangeの未定乗数法に基づく厳密な説明もできるが、ここではより直感的かつ実践的なアプローチで説明しよう。

まず結論から言ってしまうと、$\boldsymbol u_i$ はデータに関する分散共分散行列

\begin{align} S = \frac{1}{N} \sum_n (\boldsymbol x_n - \overline{\boldsymbol x})(\boldsymbol x_n - \overline{\boldsymbol x})^\top \end{align}

の固有ベクトルを、固有値が大きい順に決めればよい。

\begin{align} S\boldsymbol u_i = \lambda_i \boldsymbol u_i,\quad \lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_D \geq 0 \end{align}

ただし $\overline{\boldsymbol x} = \frac{1}{N} \sum_n \boldsymbol x_n$ である。なぜこのような方法でうまくいくのか、理由を説明しよう。まず $\boldsymbol u_i$ を $S$ の固有ベクトルに取ることで、 $S$ は $\boldsymbol u_i$ の基底から見ると対角化できる。

\begin{align} \boldsymbol u_i^\top S \boldsymbol u_j = \lambda_i \delta_{ij} \label{eq:diagonalization} \end{align}

$S$ は実対称行列なので、異なる固有値の固有ベクトルは自動的に直交するように選ばれている。もし $\lambda_i = \lambda_j$ となるような $i \neq j$ がある場合でも、$\boldsymbol u_i$ と $\boldsymbol u_j$ が直交するように基底を選べば問題はない。こうして $S$ の固有ベクトル $\boldsymbol u_i$ は正規直交基底にでき、$(\ref{eq:diagonalization})$ のように書ける。

また、$S$ は半正定値行列なので、すべての固有値は非負である。よって、$S$ の固有ベクトル $\boldsymbol u_i$ は、固有値の大きい順に並べることができる。

\begin{align} \lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_D \geq 0 \end{align}

$S$ は分散共分散行列だったので、対角化された $S$ を見ると、$\boldsymbol u_i$ の方向に沿った分散がまさに $\lambda_i$ に他ならない。したがって、$\boldsymbol u_1$ はデータの分散が最も大きい方向を表すベクトルである。$\boldsymbol u_2$ は $\boldsymbol u_1$ と直交する方向の中で、分散が最も大きい方向を表すベクトルである。以下同様に、$\boldsymbol u_M$ は $\boldsymbol u_1,\cdots,\boldsymbol u_{M-1}$ と直交する方向の中で、分散が最も大きい方向を表すベクトルである。

以下、$\boldsymbol u_1, \ldots, \boldsymbol u_M$ を主成分方向と呼ぶ。また $\boldsymbol u_i^\top \boldsymbol x_n$ を第 $i$ 主成分、または第 $i$ 主成分スコアと呼ぶ。$\boldsymbol x_n$ を $M$ 次元の主成分分析で情報を落とした近似値を $\tilde{\boldsymbol x}_n$ と書く。つまり

\begin{align} \tilde{\boldsymbol x}_n \simeq \sum_{i=1}^M \boldsymbol u_i (\boldsymbol u_i^\top \boldsymbol x_n) + C \end{align}

である。

第 $M$ 主成分までを使うとき、その $M$ 個の主成分が全体の分散のうちどれだけを説明しているかを表す量を累積寄与率という。第 $i$ 主成分だけの割合は寄与率と呼び、

\begin{align} \style{font-family:inherit;}{\text{第 }} i \style{font-family:inherit;}{\text{ 主成分の寄与率}} &= \frac{\lambda_i}{\lambda_1+\lambda_2+\cdots+\lambda_D},\\ \style{font-family:inherit;}{\text{第 }} M \style{font-family:inherit;}{\text{ 主成分までの累積寄与率}} &= \frac{\lambda_1+\cdots+\lambda_M}{\lambda_1+\lambda_2+\cdots+\lambda_D} \end{align}

で定義する。寄与率が大きい主成分ほど、元のデータのばらつきをよく説明している。したがって、低次元に落とすときには、累積寄与率を見ながら「どこまで主成分を残すか」を判断することが多い。

主成分分析の具体例

まずは $2$ 次元正規分布で考える。平均は簡単のため $\boldsymbol 0$ とし、分散共分散行列だけを見る。

例1：もともと軸がそろっている場合

まず、確率変数 $\boldsymbol x$ が

\begin{align} \boldsymbol x &\sim \mathcal{N}\left( \begin{pmatrix} 0\\ 0 \end{pmatrix}, \begin{pmatrix} 100 & 0\\ 0 & 1 \end{pmatrix} \right) \end{align}

に従うとする。この正規分布の分散共分散行列はすでに対角行列なので、固有値と固有ベクトルはすぐに分かる。

\begin{align} \lambda_1=100,\quad \boldsymbol u_1=\begin{pmatrix}1\\0\end{pmatrix}, \qquad \lambda_2=1,\quad \boldsymbol u_2=\begin{pmatrix}0\\1\end{pmatrix} \end{align}

つまり第 $1$ 主成分は $x$ 軸方向、第 $2$ 主成分は $y$ 軸方向である。第 $1$ 主成分の寄与率は

\begin{align} \frac{\lambda_1}{\lambda_1+\lambda_2} = \frac{100}{101} \simeq 0.990 \end{align}

なので、$x$ 軸方向だけで分散の約 $99.0\%$ を説明できる。この例はほとんど自明である。

ただし、ここで注意すべきなのは、単位やスケールが違う変数をそのまま使うと、分散の大きい変数が主成分を支配してしまうことである。たとえば片方が「円」、もう片方が「割合」のような場合、分散の大きさは情報量の大きさではなく単位の違いを反映しているだけかもしれない。そのため、単位の異なる変数を同列に扱う場合には、前処理としてデータを標準化し、相関行列を見ることが重要である。

例2：強く相関している場合

次に、確率変数 $\boldsymbol x$ が

\begin{align} \boldsymbol x &\sim \mathcal{N}\left( \begin{pmatrix} 0\\ 0 \end{pmatrix}, \begin{pmatrix} 1 & 0.95\\ 0.95 & 1 \end{pmatrix} \right) \end{align}

に従う場合を考える。$2$ つの変数の分散はどちらも $1$ だが、相関がかなり強い。この正規分布の分散共分散行列について、

\begin{align} \begin{pmatrix} 1 & 0.95\\ 0.95 & 1 \end{pmatrix} \begin{pmatrix} 1\\ 1 \end{pmatrix} &= 1.95 \begin{pmatrix} 1\\ 1 \end{pmatrix},\\ \begin{pmatrix} 1 & 0.95\\ 0.95 & 1 \end{pmatrix} \begin{pmatrix} 1\\ -1 \end{pmatrix} &= 0.05 \begin{pmatrix} 1\\ -1 \end{pmatrix} \end{align}

より、正規化した固有ベクトルは

\begin{align} \boldsymbol u_1 = \frac{1}{\sqrt{2}} \begin{pmatrix} 1\\ 1 \end{pmatrix}, \qquad \boldsymbol u_2 = \frac{1}{\sqrt{2}} \begin{pmatrix} 1\\ -1 \end{pmatrix} \end{align}

である。第 $1$ 主成分は「$2$ つの変数が一緒に増える方向」、第 $2$ 主成分は「片方が増えて片方が減る方向」を表している。寄与率を計算すると、

\begin{align} \style{font-family:inherit;}{\text{第 }}1\style{font-family:inherit;}{\text{ 主成分の寄与率}} &= \frac{1.95}{1.95+0.05} = 0.975,\\ \style{font-family:inherit;}{\text{第 }}2\style{font-family:inherit;}{\text{ 主成分の寄与率}} &= \frac{0.05}{1.95+0.05} = 0.025 \end{align}

となる。つまり、$2$ つの変数は見かけ上は $2$ 次元だが、ばらつきのほとんどは $\boldsymbol u_1$ 方向に乗っており、$1$ 次元で約 $97.5\%$ を説明できる。幾何的には、点の雲が右上がりの細長い楕円として分布している状況である。このとき主成分分析は、その楕円の長軸方向を第 $1$ 主成分として見つけている。

例3：見かけより低次元なデータ

恐らくもっとも教育的なのは、見かけ上は $D$ 次元だが、本当は $M$ 次元の構造に誤差が乗っているだけ、という状況だろう。例えば、各データ $\boldsymbol x_n \in \mathbb{R}^D$ が

\begin{align} \boldsymbol x_n &= A\boldsymbol z_n + \boldsymbol\varepsilon_n \end{align}

で生成されているとする。ここで $\boldsymbol z_n \in \mathbb{R}^M$ は本質的な低次元の変数、$A$ は $M$ 次元の情報を $D$ 次元の観測値に埋め込む行列、$\boldsymbol\varepsilon_n$ は小さな誤差確率変数である。つまり、観測値は $D$ 個の成分を持っているが、主要な変動は $M$ 次元の部分空間の中にある。

このとき、誤差 $\boldsymbol\varepsilon_n$ が十分小さければ、分散共分散行列の固有値は最初の $M$ 個だけが大きく、残りの $D-M$ 個は小さくなる。したがって、寄与率を見ると

\begin{align} \lambda_1+\cdots+\lambda_M \gg \lambda_{M+1}+\cdots+\lambda_D \end{align}

という形が現れる。これは「見かけは $D$ 次元だが、情報の大部分は $M$ 次元で説明できる」というサインである。主成分分析は、このような低次元の部分空間をデータから見つける方法だと考えると分かりやすい。

たとえば、画像データやアンケートデータでは、観測される変数の数は多くても、実際には「明るさ」「大きさ」「好みの強さ」のような少数の主要な変動方向が全体を動かしていることがある。主成分分析では、固有値と寄与率を見て、その少数の軸をどこまで残せばよいかを判断する。

例4：相関行列で見る場合

最後に、変数の単位が違う場合を考えよう。たとえば、あるデータに「身長」と「体重」と「握力」が含まれているとする。これらはどれも体格に関係する変数なので互いに正の相関を持ちやすいが、単位も分散の大きさも違う。このまま分散共分散行列で主成分分析を行うと、分散の大きい変数の影響が強く出る。

そこで各変数を標準化し、相関行列に対して主成分分析を行う。たとえば第 $1$ 主成分が

\begin{align} z_1 &= 0.58\,\style{font-family:inherit;}{\text{身長}} + 0.57\,\style{font-family:inherit;}{\text{体重}} + 0.58\,\style{font-family:inherit;}{\text{握力}} \end{align}

のように得られたなら、これは「総合的な体格の大きさ」を表す軸と解釈できる。一方で、第 $2$ 主成分が「体重」に正、「身長」に負の係数を持つなら、「同じ体格の中で相対的に体重が重いか軽いか」を表す軸として読めるかもしれない。主成分分析では、固有値や寄与率だけでなく、各主成分の係数を見て、その軸が何を表しているかを解釈することが大切である。

因子分析

主成分分析は、観測された変数のばらつきをよく説明する新しい軸を作る方法だった。一方で、因子分析は、観測された多くの変数の背後に、直接は観測できない少数の共通因子があると考える方法である。つまり、変数どうしが相関している理由を「背後に同じ因子が効いているから」と説明しようとする。本節では、因子分析の考え方について深堀りしない程度に説明する。

まずはデータではなく、統計モデルから考えよう。$D$ 次元の確率ベクトルを $\boldsymbol X=(X_1,\ldots,X_D)^\top$ とし、その平均を $\boldsymbol \mu=(\mu_1,\ldots,\mu_D)^\top$ とする。因子分析では、$\boldsymbol X$ が $K$ 次元の共通因子 $\boldsymbol Z=(Z_1,\ldots,Z_K)^\top$ と、各変数に固有の誤差 $\boldsymbol \varepsilon=(\varepsilon_1,\ldots,\varepsilon_D)^\top$ によって生成されると考える。

\begin{align} \boldsymbol X &= \boldsymbol \mu + \Lambda \boldsymbol Z + \boldsymbol \varepsilon \end{align}

ここで $\Lambda$ は $D \times K$ 行列であり、因子負荷量行列と呼ばれる。第 $i$ 行第 $k$ 列の成分 $\lambda_{ik}$ は、$i$ 番目の観測変数 $X_i$ に $k$ 番目の因子 $Z_k$ がどれだけ効いているかを表す量である。ただし $i=1,\ldots,D$、$k=1,\ldots,K$ である。

まずは、因子どうしが無相関である直交因子モデルとして書く。典型的には、共通因子と誤差について

\begin{align} E[\boldsymbol Z] &= \boldsymbol 0,\quad \mathrm{Var}(\boldsymbol Z)=1,\\ E[\boldsymbol \varepsilon] &= \boldsymbol 0,\quad \mathrm{Var}(\boldsymbol \varepsilon)=\Psi \end{align}

と仮定する。ただし $\Psi$ は対角行列である。$\Psi$ が対角行列であるというのは、各変数に固有の誤差どうしは相関しない、という仮定である。また、共通因子 $\boldsymbol Z$ と固有誤差 $\boldsymbol \varepsilon$ は無相関とする。

このとき、観測変数 $\boldsymbol X$ の分散共分散行列は

\begin{align} \mathrm{Var}(\boldsymbol X) &= \Lambda \Lambda^\top + \Psi \end{align}

となる。つまり、分散共分散行列を「共通因子で説明できる部分」$\Lambda \Lambda^\top$ と、「各変数に固有の部分」$\Psi$ に分けて考えるのである。

ここで、$i$ 番目の変数の分散だけを見ると

\begin{align} \mathrm{Var}(X_i) &= \lambda_{i1}^2+\lambda_{i2}^2+\cdots+\lambda_{iK}^2+\psi_i \end{align}

と書ける。右辺の $\lambda_{i1}^2+\cdots+\lambda_{iK}^2$ は、共通因子によって説明される分散であり、共通性と呼ばれる。一方、$\psi_i$ はその変数だけに固有のばらつきであり、独自性と呼ばれる。

ここまでが因子分析の統計モデルである。実際には、$\Lambda,\ \Psi$ は未知であり、$N$ 個の観測データ $(\boldsymbol x_n)_{n=1,2,\cdots,N}$ からこれらを推定する。これまでと同様に、データ（標本）の平均、分散共分散行列を $\overline{\boldsymbol x},\ S$ と書く。そして、モデル上の分散共分散行列 $\Lambda \Lambda^\top+\Psi$ が、観測データから得られた $S$ にできるだけ近くなるように、$\Lambda$ と $\Psi$ を推定する。

\begin{align} S \simeq \Lambda \Lambda^\top+\Psi \end{align}

推定方法はいくつかある。代表的には、主因子法のように共通性を仮に置いて固有値分解を繰り返す方法や、$\boldsymbol X$ に正規分布を仮定して尤度を最大にする最尤法が用いられる。いずれの場合も、観測された $S$ とモデル $\Lambda \Lambda^\top+\Psi$ のズレが小さくなるように、因子負荷量 $\Lambda$ と独自分散 $\Psi$ を調整していく、と考えればよい。

また、因子数 $K$ は自由に大きくできるわけではない。標本分散共分散行列 $S$ には、対称性を考えると独立な成分が

\begin{align} \frac{D(D+1)}{2} \end{align}

個ある。一方で、$\Lambda$ には $DK$ 個、対角行列 $\Psi$ には $D$ 個のパラメータがある。ただし因子軸の回転によって同じ分散共分散構造を表せるので、その分の自由度 $K(K-1)/2$ は差し引いて考える。したがって、モデルの自由度は

\begin{align} \frac{D(D+1)}{2} - \left\{ DK + D - \frac{K(K-1)}{2} \right\} \end{align}

である。この値が小さすぎたり負になったりする場合、データの分散共分散構造に対して因子数 $K$ が大きすぎる可能性がある。

主成分分析との違いを整理しておこう。主成分分析では、観測変数の分散をなるべくよく説明する軸を作る。したがって、主成分は観測変数の線形結合として定義される。一方、因子分析では、観測変数の背後にある潜在変数を仮定し、その潜在変数が観測変数を生み出していると考える。

\begin{align} \style{font-family:inherit;}{\text{主成分分析}} &:\quad \boldsymbol z_n = U^\top \boldsymbol x_n,\\ \style{font-family:inherit;}{\text{因子分析}} &:\quad \boldsymbol X = \boldsymbol \mu + \Lambda \boldsymbol Z + \boldsymbol \varepsilon \end{align}

前者は「観測変数から新しい軸を作る」方法であり、後者は「潜在因子から観測変数が生成される」と見る方法である。この違いはかなり重要である。主成分分析では分散の大きさそのものをよく説明することが目的だが、因子分析では変数間の相関構造を、少数の共通因子で説明することが目的である。

まとめると、因子分析は「たくさんの観測変数の相関を、少数の潜在因子で説明する」方法である。主成分分析が低次元にデータを眺め直すための幾何的な方法だとすれば、因子分析は観測変数の背後にある構造をモデルとして説明するための方法である。

因子分析では、推定された因子負荷量をそのまま読むよりも、因子軸を回転したほうが解釈しやすいことが多い。まず直交回転を考えよう。$Q$ を直交行列とすると、$\Lambda$ を $\Lambda Q$ に置き換えても

\begin{align} (\Lambda Q)(\Lambda Q)^\top &= \Lambda QQ^\top \Lambda^\top = \Lambda \Lambda^\top \end{align}

なので、共通因子で説明される分散共分散構造は変わらない。つまり、因子負荷量は一意には決まらず、同じモデルを別の因子軸で表せる。そこで、因子負荷量が解釈しやすくなるように因子軸を回転する。

回転には、大きく分けて直交回転と斜交回転がある。

直交回転：因子どうしが無相関のままになるように回転する方法である。因子軸は直角に交わるので、各因子を独立した評価軸として解釈しやすい。代表例はバリマックス回転で、各変数がなるべく少数の因子に強く反応するようにして、因子負荷量の表を読みやすくする。
斜交回転：因子どうしの相関を許して回転する方法である。因子軸は必ずしも直角に交わらないが、現実の心理的・社会的な因子は互いに相関していることも多いので、より自然な解釈になることがある。斜交回転では、単に $\Lambda\Lambda^\top$ を保つというよりも、因子間の相関行列を含めて同じ共分散構造を表す、と考える。代表例はプロマックス回転で、たとえば「価格感度」と「品質志向」のように、少し相関していそうな因子を扱いやすい。

どちらを使うべきかは、因子どうしを独立した軸として見たいのか、それとも因子間の相関も認めて自然な構造を見たいのかによって変わる。

例：マーケティングの因子分析

マーケティングでは、顧客アンケートの背後にある潜在的な評価軸を探すために因子分析が使われる。たとえば、アンケート結果を用いて、観光客が地域の魅力・経験価値をどのように評価しているかを因子分析してみよう。なお、本節の説明は下記出典の論文を参考にしているが、ここでは発想を単純化して簡単に説明する。

出典：

那須一貴「地域の魅力度と地域での経験価値が観光客の満足度と再来訪意向に及ぼす影響に関する考察」文教大学国際学部紀要第28巻第2号, 2018年。PDFはこちら。

たとえば、観光客に「楽しい気持ちになれそうか？」「他の地域にはない経験が得られそうか」といったアンケートに答えてもらったとしよう。その回答を因子分析すると、たくさんの質問項目の背後にある少数の評価軸を取り出せる。ここでは、地域の経験価値に関する $5$ つの質問項目から、「ワクワク感」と「特別感」の $2$ 因子だけを抜き出した表として見てみよう。つまり $D=5$、$K=2$ である。

\begin{align} \begin{array}{c|cc} & \style{font-family:inherit;}{\text{ワクワク感}} & \style{font-family:inherit;}{\text{特別感}}\\ \hline \style{font-family:inherit;}{\text{時代のトレンドや話題性が感じられる}} & 0.928 & 0.002\\ \style{font-family:inherit;}{\text{なじみがあり、安心できる}} & 0.800 & -0.170\\ \style{font-family:inherit;}{\text{楽しい気持ちになれる}} & 0.739 & 0.197\\ \style{font-family:inherit;}{\text{他の地域にはない経験が得られる}} & -0.087 & 1.033\\ \style{font-family:inherit;}{\text{自分を見つめなおすことができる}} & -0.230 & 0.806 \end{array} \end{align}

「時代のトレンドや話題性が感じられる」「なじみがあり、安心できる」「楽しい気持ちになれる」は「ワクワク感」に強く効いている。一方、「他の地域にはない経験が得られる」「自分を見つめなおすことができる」は「特別感」に強く効いている。このように、複数のアンケート項目を少数の因子名でまとめて解釈できる。

さらに、各観光地の因子スコアを推定すれば、観光地ごとの位置づけを次のように図示できる。因子スコアとは、ある観光地 $\boldsymbol x_n$ が各因子をどれくらい持っているかを表す推定値である。代表的な回帰法では、推定済みの $\hat{\Lambda}$ と $\hat{\Psi}$ を用いて

\begin{align} \hat{\boldsymbol z}_n &= \hat{\Lambda}^\top (\hat{\Lambda}\hat{\Lambda}^\top+\hat{\Psi})^{-1} (\boldsymbol x_n-\overline{\boldsymbol x}) \end{align}

のように計算する。今回の例では $\hat{\boldsymbol z}_n$ の第 $1$ 成分を「ワクワク感」、第 $2$ 成分を「特別感」と読んでいる。ただし、これは因子スコアの代表的な推定法の一つである。因子スコアの推定法には Bartlett 法などもあり、方法によって値は少し変わる。

横軸を「ワクワク感」、縦軸を「特別感」とすると、各観光地がどの経験価値に訴えているかを比較しやすい。なお、次の図は（元論文とは関係なく）因子スコアの使い方を示す模式図である。

たとえば「話題の温泉街A」はワクワク感が強く、特別感は低めである。一方、「内省の宿B」はワクワク感よりも特別感で評価されている。「非日常体験D」は両方が高く、楽しい気持ちと他地域にはない経験の両方を持っている観光地だと読める。このように、因子分析はアンケート項目を減らすだけでなく、観光地を解釈しやすい軸の上に配置するためにも使える。

ただし、この解釈は機械的に決まるものではない。因子分析では、数学的に推定された因子負荷量を見て、人間がその因子に意味を与える必要がある。また、因子の符号は本質的ではない。ある因子の負荷量をすべて $-1$ 倍しても、モデルとしては同じ意味を持つ。

品質管理

品質管理では、製品や工程が安定しているか、また規格を満たす能力があるかを統計的に確認する。ここでは代表的な道具として、管理図、工程能力指数、不良率を見ていこう。

まず、工程が安定しているかどうかを見るために使うのが管理図である。たとえば、ある製品の寸法を一定時間ごとに測定し、その平均値を時系列でプロットする。もし工程が安定していれば、測定値は平均のまわりでランダムにばらつくはずである。一方、突然大きく外れたり、同じ方向に偏り続けたりする場合は、工程に何らかの異常が起きている可能性がある。

代表的な管理図では、中心線を工程平均 $\mu$ とし、その上下に管理限界を引く。安定状態における標準偏差を $\sigma$ とすると、よく用いられるのは $3\sigma$ 管理限界である。

\begin{align} \mathrm{UCL} &= \mu + 3\sigma,\\ \mathrm{CL} &= \mu,\\ \mathrm{LCL} &= \mu - 3\sigma \end{align}

ここで $\mathrm{UCL}$ は上方管理限界（Upper Control Limit）、$\mathrm{LCL}$ は下方管理限界（Lower Control Limit）、$\mathrm{CL}$ は中心線（Center Line）である。観測値が管理限界の外に出たときは、偶然のばらつきだけでは説明しにくい異常として扱う。ただし、管理図は「規格に合っているか」を見る図ではなく、「工程が安定しているか」を見る図である。この違いは重要である。

次に、工程が規格に対してどれだけ余裕を持っているかを表す量として、工程能力指数がある。製品にはしばしば上限規格 $\mathrm{USL}$ と下限規格 $\mathrm{LSL}$ が決められている。たとえば、ある寸法が $10.00 \pm 0.05$ mm でなければならないなら、$\mathrm{USL}=10.05$、$\mathrm{LSL}=9.95$ である。

工程のばらつきが標準偏差 $\sigma$ で表されるとき、工程能力指数 $C_p$ は

\begin{align} C_p &= \frac{\mathrm{USL}-\mathrm{LSL}}{6\sigma} \end{align}

で定義される。分母の $6\sigma$ は、平均の上下 $3\sigma$ を合わせた幅である。つまり $C_p$ は、規格幅が工程の自然なばらつき幅に対してどれだけ広いかを表している。$C_p$ が大きいほど、ばらつきに対して規格幅に余裕がある。

ただし、$C_p$ は工程平均が規格の中心にあることを暗黙に期待している。実際には、ばらつきが小さくても平均が規格の端に寄っていれば、不良品は出やすくなる。そこで、平均のずれも考慮した工程能力指数として $C_{pk}$ を用いる。

\begin{align} C_{pk} &= \min\left( \frac{\mathrm{USL}-\mu}{3\sigma}, \frac{\mu-\mathrm{LSL}}{3\sigma} \right) \end{align}

$C_{pk}$ は、平均 $\mu$ から見て、上限規格と下限規格のどちらに近いかを考慮する。工程平均が規格の中心からずれると、片側の余裕が小さくなるので、$C_{pk}$ は $C_p$ より小さくなる。したがって、$C_p$ は「ばらつきの小ささ」、$C_{pk}$ は「ばらつきの小ささと中心のずれ」を見る指標だと考えると分かりやすい。

最後に、不良率を考えよう。工程から出てくる製品の特性値 $X$ が正規分布

\begin{align} X &\sim \mathcal{N}(\mu,\sigma^2) \end{align}

に従うと仮定する。このとき、不良品とは $X$ が規格外に出る製品である。したがって、不良率は

\begin{align} P(\style{font-family:inherit;}{\text{不良}}) &= P(X < \mathrm{LSL}) + P(X > \mathrm{USL})\\ &= \Phi\left(\frac{\mathrm{LSL}-\mu}{\sigma}\right) + 1-\Phi\left(\frac{\mathrm{USL}-\mu}{\sigma}\right) \end{align}

で計算できる。ただし $\Phi$ は標準正規分布の累積分布関数である。この式から分かるように、不良率を下げるには、$\sigma$ を小さくしてばらつきを抑えることと、$\mu$ を規格の中心に近づけることが重要である。

たとえば、規格が $\mu \pm 3\sigma$ に相当する幅であれば、正規分布では規格外に出る確率は約 $0.27\%$ である。一方、規格が $\mu \pm 6\sigma$ に相当する幅であれば、規格外に出る確率はきわめて小さくなる。このように、管理図は工程の安定性を監視し、工程能力指数は規格に対する余裕を測り、不良率は実際に規格外へ出る確率を評価するために使われる。

ハザード関数

ハザード関数は、生存分析や信頼性工学で使われる概念である。ある対象が「生きている」状態から「死ぬ」状態に移る確率を時間とともに表す関数である。たとえば、機械部品の故障率や患者の死亡率などをモデル化するために用いられる。

ハザード関数の定義は難しくない。簡単に言うと、ハザード関数 $\lambda(t)$ は、単位時間あたりに製品が故障する確率である。より正確には、時刻 $t$ で製品が故障していないとき、次の瞬間 $\mathrm{d}t$ の間に故障する確率を $\lambda(t)\mathrm{d}t$ とする。なお、ハザード関数は「瞬間的に死ぬ確率」と説明されることがあるが、これは正確ではない。より正確には、$\lambda(t)\mathrm{d}t$ の形で理解されるものである。

さて、時刻 $t=0$ で製品が故障しない状態で使われ始めたとしよう。このとき、時刻 $t$ に故障している確率を $F(t)$ と書く。すると、時刻 $t$ で製品が故障していない確率は $1 - F(t)$ である。これに、次の瞬間 $\mathrm{d}t$ の間に故障する確率 $\lambda(t)\mathrm{d}t$ をかけて、 $\mathrm{d}t$ の間に初めて故障する確率 $\mathrm{d}F(t)$ は

\begin{align} \mathrm{d}F(t) &= \big( 1 - F(t) \big) \lambda(t) \mathrm{d}t \label{eq:hazard-function-differential} \end{align}

と書ける。あとはこれを微分方程式としてみて解くだけである。両辺を $1 - F(t)$ で割って積分することにより

\begin{align} \log(1 - F(t)) &= \int_0^t \lambda(t') \, \mathrm{d}t' \end{align}

つまり

\begin{align} F(t) &= 1 - \exp\left(-\int_0^t \lambda(t') \, \mathrm{d}t'\right) \end{align}

が得られる。確率密度関数 $f(t) = \frac{\mathrm{d}F(t)}{\mathrm{d}t}$ を求めると、

\begin{align} f(t) &= \lambda(t) \exp\left(-\int_0^t \lambda(t') \, \mathrm{d}t'\right) \end{align}

である。つまり、ハザード関数 $\lambda(t)$ を知れば、故障する確率 $F(t)$ や確率密度関数 $f(t)$ を計算できる。逆に、$F(t)$ や $f(t)$ を知れば、$\lambda(t)$ を求めることもできる。これは $(\ref{eq:hazard-function-differential})$ により

\begin{align} \lambda(t) &= \frac{f(t)}{1 - F(t)} \end{align}

と書けるからである。ハザード関数は、時間とともに故障率がどのように変化するかを表す重要な指標である。たとえば、製品の初期不良が多い場合は、$\lambda(t)$ は最初に高くなり、その後下がる形になることが多い。一方で、摩耗による故障が多い場合は、$\lambda(t)$ は時間とともに増加する形になることが多い。$\lambda(t)$ が時刻 $t$ によらない一定値 $\lambda$ である場合は指数分布に等しく、原子核の崩壊などがこれに従う。半減期は

\begin{align} T_{1/2} &= \frac{1}{\lambda}\log 2 \end{align}

である。分母に $\lambda$ があることから、ハザード関数が大きいほど半減期は短くなることが分かる。

一方で、$\lambda(t)$ が

\begin{align} ax^b = \int_0^t \lambda(t') \, \mathrm{d}t' \end{align}

のような形で表されることもある。これはワイブル分布と呼ばれる分布に対応し、確率密度関数は

\begin{align} f(t) &= ab t^{b-1} \exp(-a t^b) \end{align}

となる。ワイブル分布は、$\lambda(t)$ が時間とともに増加したり減少したりする場合に適している。たとえば、製品の初期不良が多い場合は $b<1 $、摩耗による故障が多い場合は $b >1$ となることが多い。

ポアソン過程

窓口に顧客がつねに一定の確率でランダムにやってくるとき、時刻 $t$ までの間にやってくる顧客数を $N_t$ とする。このように、時間とともに進む確率変数は一般に確率過程と呼ばれ、特にこのような設定の確率過程はポアソン過程と呼ばれる。ここでは、確率過程について数学的に厳密な取り扱いは行わないし、そのつもりもない。その上で、このような分布について考えよう。

ポアソン過程の定義はこうである。$\lambda$ を正の定数とする。$N_t$ は $t \geq 0$ に対して定義された確率過程であって、次の条件を満たすとき、$N_t$ はパラメータ $\lambda$ のポアソン過程であると言う。

$N_0 = 0$ である。
確率過程 $N_t$ は独立増分である。つまり、任意の $0 \leq t_1 < t_2 < \cdots < t_n$ に対して、$N_{t_1}-N_0, N_{t_2}-N_{t_1}, \cdots, N_{t_n}-N_{t_{n-1}}$ は互いに独立である。
確率過程 $N_t$ は定常増分である。つまり、$N_{t+s}-N_s$ の分布は $s$ によらず、時間の長さ $t$ だけで決まる。
微小時間 $\mathrm{d}t$ の間に顧客は $1$ 人だけやってくるとみなしてよく、その確率はほぼ $\lambda \mathrm{d}t $ である。より正確には、$P(N_{t+\mathrm{d}t} - N_t = 1) = \lambda \mathrm{d}t + o(\mathrm{d}t)$ 、$P(N_{t+\mathrm{d}t} - N_t \geq 2) = o(\mathrm{d}t)$ である。

それでは $P(N_t = k)$ の値を求めよう。いくつか説明方法はあるが、最も手っ取り早いのは $[0,t]$ を $n$ 個の微小時間区間 $t = n \Delta t$ に分けて、二項分布で近似してしまうことである。$n$ を大きく取れば、各区間 $\Delta t$ にやってくる顧客の人数は高々 $1$ 人であり、その確率が $\lambda \Delta t$ となる。よって、

\begin{align} P(N_t = k) &= \lim_{n \to \infty} {}_n \mathrm{C}_k (\lambda \Delta t)^k (1 - \lambda \Delta t)^{n-k} \end{align}

である。$\Delta t = t/n$ を代入して、$n \to \infty$ （$\Delta t \to 0$）の極限を取ると、以前計算したように

\begin{align} P(N_t = k) &= \frac{(\lambda t)^k}{k!} \exp(-\lambda t) \end{align}

が得られる。つまり、$N_t$ はパラメータ $\lambda t$ のポアソン分布に従うことが分かる。したがって、ポアソン過程では、時刻 $t$ までの間にやってくる顧客数は平均 $\lambda t$ のポアソン分布に従うことになる。

ここまでで「時刻 $t$ までに何回起こるか」が分かった。ポアソン過程のおもしろいところは、ここで終わらないところである。同じ考え方で、任意の時間区間 $(s,t]$ だけを切り出しても、そこに起こる回数は区間の長さ $t-s$ だけで決まる。

\begin{align} N_t - N_s &\sim \mathrm{Po}(\lambda (t-s)) \end{align}

さらに、互いに重ならない区間で起こる回数は独立である。つまり、ポアソン過程では「どの時刻から数え始めたか」ではなく、「どれだけの長さを観測したか」が本質である。この見方に慣れると、待ち時間の分布も自然に出てくる。

次に、人数ではなく「何時に来るか」を考える。最初の到着時刻を $T_1$ とする。$T_1 > t$ であることは、時刻 $t$ までに誰も来ないことと同じなので、

\begin{align} P(T_1 > t) &= P(N_t = 0) = \exp(-\lambda t) \end{align}

である。したがって $T_1$ はパラメータ $\lambda$ の指数分布に従う。ここが、ポアソン過程を理解する上でかなり大事である。「回数」を見ればポアソン分布、「次に起こるまでの時間」を見れば指数分布になる。

第 $n$ 番目の到着時刻を $T_n$ とし、到着間隔を $S_n = T_n - T_{n-1}$（ただし $T_0=0$）と書くと、到着間隔は互いに独立な指数分布になる。

\begin{gather} S_1,S_2,\ldots \overset{\mathrm{i.i.d.}}{\sim} \mathrm{Exp}(\lambda), \\ T_n = S_1+\cdots+S_n \sim \mathrm{Gamma}(n,\lambda) \end{gather}

ここで $\mathrm{Gamma}(n,\lambda)$ は、形状母数 $n$、率母数 $\lambda$ のガンマ分布である。第 $n$ 番目の到着まで待つということは、指数分布に従う待ち時間を $n$ 個足すということなので、ガンマ分布が出てくる。密度は

\begin{align} f_{T_n}(t) &= \frac{\lambda^n t^{n-1}}{(n-1)!}\exp(-\lambda t) \qquad (t>0) \end{align}

である。また、$T_n \leq t$ は「時刻 $t$ までに少なくとも $n$ 回起こる」ことと同じなので、

\begin{align} P(T_n \leq t) &= P(N_t \geq n) \end{align}

である。この対応は、ポアソン分布とガンマ分布をつなぐ重要な関係である。問題を解くときも、「回数で考えるべきか、待ち時間で考えるべきか」を切り替えるだけで見通しがよくなることが多い。

条件付き分布もよく問われる。時刻 $t$ までにちょうど $n$ 回起こったことが分かっているとする。このとき、その $n$ 個の到着時刻は、区間 $(0,t)$ の中に一様にばらまかれた点を小さい順に並べたものと同じ分布をもつ。これは少し意外に見えるが、$N_t=n$ と分かった後では、どの場所に起こりやすいという偏りは残らないのである。$U_1,\ldots,U_n$ を互いに独立な一様分布 $\mathrm{Unif}(0,t)$ に従う確率変数とし、その順序統計量を $U_{(1)} < \cdots < U_{(n)}$ と書くと、

\begin{align} (T_1,\ldots,T_n)\mid (N_t=n) &\overset{d}{=} (U_{(1)},\ldots,U_{(n)}) \end{align}

である。ここで $\overset{d}{=}$ は、確率変数そのものが同じという意味ではなく、同じ分布をもつという意味である。つまり、時刻 $t$ までに $n$ 回起こったと分かったときの到着時刻の並びは、$n$ 個の一様乱数を小さい順に並べたものと同じ分布をもつ。ここから、たとえば $0< s < t$ に対して、時刻 $s$ までに起こった回数の条件付き分布は

\begin{align} N_s \mid (N_t=n) &\sim \mathrm{Bin}\left(n,\frac{s}{t}\right) \end{align}

となる。$t$ までに起こった $n$ 回のうち、それぞれが区間 $(0,s]$ に入る確率が $s/t$ である、と見れば二項分布になるのは自然である。また、第 $k$ 到着時刻については

\begin{align} \frac{T_k}{t}\mid (N_t=n) &\sim \mathrm{Beta}(k,n-k+1), \\ E[T_k\mid N_t=n] &= \frac{k}{n+1}t \end{align}

が成り立つ。

ポアソン過程には、分割と重ね合わせという扱いやすい性質もある。たとえば、来店客を独立に確率 $p$ で「会員」、確率 $1-p$ で「非会員」に分類するとしよう。このとき、会員の到着回数 $N_t^{(A)}$ と非会員の到着回数 $N_t^{(B)}$ は、それぞれ

\begin{align} N_t^{(A)} &\sim \mathrm{Po}(\lambda p t), \\ N_t^{(B)} &\sim \mathrm{Po}(\lambda (1-p)t) \end{align}

に従い、しかも互いに独立である。これはポアソン過程の分割と呼ばれる。もとの到着をランダムに色分けしただけなのに、それぞれがまたポアソン過程になるというのが気持ちのよい性質である。

逆に、独立なポアソン過程 $N_t^{(1)},\ldots,N_t^{(m)}$ があり、それぞれの強度が $\lambda_1,\ldots,\lambda_m$ であるなら、その和

\begin{align} N_t &= N_t^{(1)}+\cdots+N_t^{(m)} \end{align}

は、強度 $\lambda_1+\cdots+\lambda_m$ のポアソン過程になる。これはポアソン過程の重ね合わせである。別々の入口から人が入ってくる状況を合計して、全体の到着過程として眺める、というイメージである。

最後に、推定の形も確認しておこう。時間区間 $[0,T]$ でポアソン過程を観測し、その間の到着回数が $N_T=n$ であったとする。このとき、$\lambda$ に関する尤度は、定数倍を無視すれば

\begin{align} L(\lambda) &\propto \exp(-\lambda T)\lambda^n \end{align}

である。したがって対数尤度は

\begin{align} \ell(\lambda) &= -\lambda T + n\log\lambda + \mathrm{const.} \end{align}

となり、最尤推定量は

\begin{align} \hat{\lambda} &= \frac{N_T}{T} \end{align}

である。これは「単位時間あたりに何回起こったか」という、直感通りの推定量になっている。均質なポアソン過程では、到着時刻そのものを観測していても、$\lambda$ について本質的な情報は総到着回数 $N_T$ に集約される。条件付き分布 $(T_1,\ldots,T_n)\mid(N_T=n)$ は $\lambda$ を含まないからである。細かい時刻のデータを持っているのに、推定では結局 $N_T/T$ に落ちる、というのは少し驚くが、均質なポアソン過程らしいきれいな結論である。

統計検定一級の問題として見るなら、$N_t$ のポアソン分布だけでなく、$T_n$ のガンマ分布、$N_s\mid N_t$ の二項分布、到着時刻の順序統計量、分割と重ね合わせ、そして $\hat{\lambda}=N_T/T$ までを一つのまとまりとして押さえておくとよい。