統計学の基礎をすぐに思い出せるような、まとめチートシートを作成しました。
「確率分布編」「推定・検定編」「線形回帰編」「その他の話題（理工学分野）」の４部構成です。本ページは「その他の話題（理工学分野）」です。
注意
あくまで個人用のメモであり、誤植や不正確な情報が含まれている可能性があります。適宜更新していきます。
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$ を用いて、なるべく近似したいとしよう。つまり、
となるような基底 $\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$ 成分の分散は、ほとんどゼロに近いであろう。
なお、実際は $\boldsymbol x_n$ の平均を引いて中心化してから主成分分析を行うことが多い。すなわち、$\boldsymbol x_n$ を $\boldsymbol x_n - \overline{\boldsymbol x}$ に置き換えて解析する。この場合、上の式は
と書き換えられる。先ほどの定数 $C$ が今回は消えていることに注意しよう。前もってデータを中心化しておくことで、主成分分析のときに定数 $C$ を考慮する必要がなくなるのである。
さて、このような基底はどのようにすれば作れるのだろうか？ これについては、誤差関数とLagrangeの未定乗数法に基づく厳密な説明もできるが、ここではより直感的かつ実践的なアプローチで説明しよう。
まず結論から言ってしまうと、$\boldsymbol u_i$ はデータに関する分散共分散行列
の固有ベクトルを、固有値が大きい順に決めればよい。
ただし $\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 - \overline{\boldsymbol x})$ を第 $i$ 主成分、または第 $i$ 主成分スコアと呼ぶ。$\boldsymbol x_n$ を $M$ 次元の主成分分析で情報を落とした近似値を $\tilde{\boldsymbol x}_n$ と書く。つまり
である。
第 $M$ 主成分までを使うとき、その $M$ 個の主成分が全体の分散のうちどれだけを説明しているかを表す量を累積寄与率という。第 $i$ 主成分だけの割合は寄与率と呼び、
で定義する。寄与率が大きい主成分ほど、元のデータのばらつきをよく説明している。したがって、低次元に落とすときには、累積寄与率を見ながら「どこまで主成分を残すか」を判断することが多い。
主成分分析の具体例
まずは $2$ 次元正規分布で考える。平均は簡単のため $\boldsymbol 0$ とし、分散共分散行列だけを見る。
例1：もともと軸がそろっている場合
まず、確率変数 $\boldsymbol x$ が
に従うとする。この正規分布の分散共分散行列はすでに対角行列なので、固有値と固有ベクトルはすぐに分かる。
\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$ が
に従う場合を考える。$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$ が
で生成されているとする。ここで $\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$ によって生成されると考える。
ここで $\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$ の分散共分散行列は
となる。つまり、分散共分散行列を「共通因子で説明できる部分」$\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$ は差し引いて考える。したがって、モデルの自由度は
である。この値が小さすぎたり負になったりする場合、データの分散共分散構造に対して因子数 $K$ が大きすぎる可能性がある。
主成分分析との違いを整理しておこう。主成分分析では、観測変数の分散をなるべくよく説明する軸を作る。したがって、主成分は観測変数の線形結合として定義される。一方、因子分析では、観測変数の背後にある潜在変数を仮定し、その潜在変数が観測変数を生み出していると考える。
前者は「観測変数から新しい軸を作る」方法であり、後者は「潜在因子から観測変数が生成される」と見る方法である。この違いはかなり重要である。主成分分析では分散の大きさそのものをよく説明することが目的だが、因子分析では変数間の相関構造を、少数の共通因子で説明することが目的である。
まとめると、因子分析は「たくさんの観測変数の相関を、少数の潜在因子で説明する」方法である。主成分分析が低次元にデータを眺め直すための幾何的な方法だとすれば、因子分析は観測変数の背後にある構造をモデルとして説明するための方法である。
因子分析では、推定された因子負荷量をそのまま読むよりも、因子軸を回転したほうが解釈しやすいことが多い。まず直交回転を考えよう。$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$ 倍しても、モデルとしては同じ意味を持つ。
品質管理
品質管理では、似ているようで目的の異なる $2$ 種類の境界を扱う。ひとつは工程の変化を見つけるための管理限界、もうひとつは製品として許容できる範囲を表す規格限界である。まず、この $2$ つを切り分けて考えよう。
管理図と管理限界
管理図の目的は、工程がいつもの安定した状態にあるかを時系列で監視することである。ここでは、過去 $1$ か月のうち $m$ 日について、毎日ほぼ同じ条件のもとで $n$ 個の製品をランダムに抽出したとしよう。基準となるこの期間では工程が安定しており、各測定値は互いに独立に、同じ正規分布に従うものとする。第 $i$ 日の $j$ 個目の測定値を $X_{ij}$ と書き、その日の標本平均と範囲を
\begin{align} \overline{X}_i &= \frac{1}{n}\sum_{j=1}^n X_{ij},\\ R_i &= \max_j X_{ij}-\min_j X_{ij} \end{align}とする。さらに、$m$ 日分の標本平均の平均と、範囲の平均を
とする。工程がこの期間を通じて安定しているなら、$\overline{\overline{X}}$ を工程平均の推定値 $\hat{\mu}$ として使える。
次に、工程のばらつきの大きさを推定する。工程標準偏差は、範囲の平均 $\overline{R}$ を用いて
\begin{align} \hat{\sigma} &= \frac{\overline{R}}{d_2} \end{align}と推定する。ここで $d_2$ は、毎日測る製品の個数 $n$ によって決まる定数であり、管理図用の係数表から選ぶ。
個々の測定値ではなく日ごとの標本平均 $\overline{X}_i$ を監視する場合、その標準誤差の推定値は $\hat{\sigma}/\sqrt{n}$ である。したがって、$\overline{X}$ 管理図の中心線と $3\sigma$ 管理限界を
と定める。$\mathrm{UCL}$ は上方管理限界、$\mathrm{LCL}$ は下方管理限界、$\mathrm{CL}$ は中心線である。新たに得られた $\overline{X}_i$ が管理限界の外に出た場合は、設備、材料、作業条件などに、通常とは異なる原因が生じた可能性を疑う。また、点が管理限界内にあっても、一方向への連続した偏りなどがあれば、工程の変化を検討する必要がある。
範囲 $R_i$ は、工程のばらつきそのものが変化していないかを見るために使う。$R$ 管理図では、毎日測る製品の個数 $n$ によって決まる定数 $D_3,D_4$ を用いて
\begin{align} \mathrm{UCL}_R &= D_4\overline{R}, & \mathrm{CL}_R &= \overline{R}, & \mathrm{LCL}_R &= D_3\overline{R} \end{align}とする。$D_3,D_4$ も管理図用の係数表から選ぶ。平均の管理図だけでなく、ばらつきの管理図も併せて確認することで、工程の状態をより適切に監視できる。
規格限界と工程能力
一方、規格限界の目的は、製品が要求された性能や寸法を満たしているかを判定することである。上側規格限界を $\mathrm{USL}$、下側規格限界を $\mathrm{LSL}$ と書く。これらは顧客要求、設計条件、法令などから定められる値であり、通常は測定データから推定するものではない。
したがって、管理限界と規格限界は一致する必要がない。日ごとの標本平均が管理限界の外に出ても、測定した個々の製品が規格外とは限らない。逆に、管理図上の点が管理限界内で推移していても、工程のばらつきが規格幅より大きければ規格外品は発生する。管理図は「工程が安定しているか」を見て、規格限界は「製品が要求を満たしているか」を見るのである。
工程が管理状態にあり、その分布が概ね正規分布とみなせ、標準偏差が時間を通じて安定しているとする。このとき、工程能力指数は
で定義される。「工程能力あり」と判断する目安として $\hat C_p\geq1.33$ がよく用いられるが、必要な水準は製品や業界によって異なる。
$\hat C_p$ は規格幅と工程の自然なばらつき幅 $6\sigma$ の比であり、工程平均の位置は考慮しない。工程平均のずれも考慮するときは、
を用いる。$\hat{C}_p$ は主に工程のばらつきの小ささを、$\hat{C}_{pk}$ はばらつきの小ささと工程平均のずれを併せて評価する。
次に、不良率を考えよう。工程から得られる製品寸法を、推定した正規分布 $\mathcal{N}(\hat{\mu},\hat{\sigma}^2)$ で近似すると、推定不良率は
\begin{align} \widehat{P}(\style{font-family:inherit;}{\text{不良}}) &= \Phi\left( \frac{\mathrm{LSL}-\hat{\mu}}{\hat{\sigma}} \right) +1- \Phi\left( \frac{\mathrm{USL}-\hat{\mu}}{\hat{\sigma}} \right) \end{align}と計算できる。ただし $\Phi$ は標準正規分布の累積分布関数である。これは正規分布という仮定に基づく推定値であり、工程が不安定な場合や分布が正規分布から大きく外れる場合には、解釈に注意が必要である。
最後に、$C_p$ の点推定値と信頼区間を考える。これらはいずれも、$n$ 回の測定によって得られた標本分散の不偏推定量 $\hat \sigma_V^2$ を用いて表現できる。まず、点推定値 $\hat{C}_p$ は \eqref{eq:cp} 式に $\hat{\sigma}=\hat \sigma_V$ を代入すれば得られる。次に、 $C_p$ の信頼区間を求めよう。まず、
\begin{align} \frac{(n-1)\hat{\sigma}_V^2}{\sigma^2} \sim \chi^2_{n-1} \end{align}である。これにより、工程能力指数の母数を $C_p = \frac{\mathrm{USL}-\mathrm{LSL}}{6\sigma}$ とすると、
\begin{align} \frac{(n-1) C_p^2}{\hat C_p^2} \sim \chi^2_{n-1} \end{align}である。これより、信頼係数 $1-\alpha$ の $C_p$ の信頼区間は
\begin{align} \chi^2_{n-1,1-\alpha/2} \leq \frac{(n-1) C_p^2}{\hat C_p^2} \leq \chi^2_{n-1,\alpha/2} \end{align}を母数 $C_p$ について解いて、
から信頼区間が得られる。
ハザード関数
ハザード関数は、生存分析や信頼性工学で使われる概念である。ある対象が「生きている」状態から「死ぬ」状態に移る確率を時間とともに表す関数である。たとえば、機械部品の故障率や患者の死亡率などをモデル化するために用いられる。
ハザード関数の定義は難しくない。簡単に言うと、ハザード関数 $\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)$ は
と書ける。あとはこれを微分方程式としてみて解くだけである。両辺を $1 - F(t)$ で割って積分することにより
\begin{align} \log(1 - F(t)) &= \int_0^t \lambda(t') \, \mathrm{d}t' \end{align}つまり
が得られる。確率密度関数 $f(t) = \frac{\mathrm{d}F(t)}{\mathrm{d}t}$ を求めると、
である。つまり、ハザード関数 $\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)$ が
のような形で表されることもある。これはワイブル分布と呼ばれる分布に対応し、確率密度関数は
となる。ワイブル分布は、$\lambda(t)$ が時間とともに増加したり減少したりする場合に適している。たとえば、製品の初期不良が多い場合は $b<1 $、摩耗による故障が多い場合は $b >1$ となることが多い。
ポアソン過程
窓口に顧客がつねに一定の確率でランダムにやってくるとき、時刻 $t$ までの間にやってくる顧客数を $N_t$ とする。このように、時間とともに進む確率変数の列 $(N_t)_{t \geq 0}$ は一般に確率過程と呼ばれ、特にこのような設定の確率過程はポアソン過程と呼ばれる。ここでは、確率過程について数学的に厳密な取り扱いに深追いしない程度に、これについて考えてみよう。
ポアソン過程の定義はこうである。$\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$ となる。より正確に言えば、$\ell = 1,\dots,n$ に対して、
- $N_t$ の独立増分性より、$N_{\ell \Delta t} - N_{(\ell-1) \Delta t}$ はすべて互いに独立である
- $N_t$ の定常増分性より、$P(N_{\ell \Delta t} - N_{(\ell-1) \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) \label{eq:poisson} \end{align}が得られる。これより、$N_t$ はパラメータ $\lambda t$ のポアソン分布に従うことが分かる。つまり、
である。
以下、$P_k(t) = P(N_t = k)$ と書く。ポアソン過程の確率分布 $P_k(t)$ は、以下の微分方程式に従う。
この微分方程式は、ポアソン過程の定義から直接導くこともできる。教科書によっては、上の微分方程式を出発点としてポアソン分布の確率分布 \eqref{eq:poisson} を導いているものもある（ここでは詳細は省略する）。
ここからは、ポアソン過程の性質についてさらに見ていこう。まず、下記が成り立つことは明らかであろう。
言い換えると、任意の時間区間 $(s,t]$ だけを切り出しても、そこに起こる回数は区間の長さ $t-s$ だけで決まる。
次に、人数ではなく、$n$ 人目が来る時刻 $S_n$ の分布を求めよう。最初に注意点として、$S_n$ は（$N_t$ と異なり）連続型確率変数なので、$P(S_n = t) = 0$ であることに注意しよう。そこで、分布関数 $P(S_n \leq t)$ を求めに行く。まず、
\begin{align} P(S_1 \leq t) &= 1 - P(S_1 > t)\\ &= 1 - P(N_t = 0) = 1 - \exp(-\lambda t) \end{align}である。これより、$S_1 \sim \mathrm{Ex}(\lambda)$ である。
次に、$R_n = S_n - S_{n-1}$ とおく。つまり、$R_n$ は $n-1$ 人目が到着してから $n$ 人目が到着するまでの時間間隔である。ただし $R_1 = S_1$ とする。このときまず、
\begin{align} S_n = R_1 + R_2 + \cdots + R_n \end{align}であり、また $R_1 \sim \mathrm{Ex}(\lambda)$ である。ポアソン過程では、到着時刻 $S_{n-1}$ より後の増分は、それ以前の過程と独立であり、もとのポアソン過程と同じ確率法則に従う。この性質を用いると、
\begin{align} P(R_n \leq t \mid S_{n-1} = s) &= P(S_n \leq s + t \mid S_{n-1} = s) \\ &= P(N_{s+t} - N_s \geq 1) \\ &= 1 - P(N_{s+t} - N_s = 0) \\ &= 1 - \exp(-\lambda t) \end{align}である。よって、任意の $n$ について、$R_n$ はそれ以前の到着過程と独立に指数分布 $\mathrm{Ex}(\lambda)$ に従う。以上より、
である。これより、$\mathrm{Ex}(\lambda) = \mathrm{Ga}(1, 1/\lambda)$ であることに注意すると、$n$ 回の和により
である。
マルコフ連鎖
有限（または可算）集合 $S$ を状態空間と呼び、その要素 $i \in S$ を状態と呼ぶ。離散時刻 $n=0,1,2,\cdots$ において、状態空間 $S$ 上の値をとる確率変数を $X_n$ とする。
ここで、次の状態 $X_{n+1}$ が、それ以前の履歴には依存せず、直前の状態 $X_n$ のみに依存する確率過程を考えよう。すなわち、
を満たす $ (X_n)_{n \geq 0} $ を考えよう。このような性質はマルコフ性と呼ばれ、上の値は状態 $i$ から $j$ への遷移確率と呼ばれる。
特に、 $i$ から $j$ への遷移確率が時刻 $n$ に依存しない場合、$p_{ij}$ を $n$ によらない定数として
と定義し、$P=(p_{ij})_{i,j \in S}$ を遷移確率行列として定義できる。上記 \eqref{eq:markov_property} および \eqref{eq:transition_probability} 式をともに満たす確率過程 $(X_n)_{n \geq 0}$ をマルコフ連鎖と呼ぶ。
時刻 $n$ における状態分布を、行ベクトルで以下のように書くことにしよう。
\begin{align} \boldsymbol{\pi}^{(n)} = \begin{pmatrix} \pi_1^{(n)} & \pi_2^{(n)} & \cdots \end{pmatrix} \end{align}ここで $\pi_i^{(n)} = P(X_n = i)$ である。これを用いると、
\begin{align} \pi_j^{(n+1)} &= P(X_{n+1} = j) \\ &= \sum_{i \in S} P(X_{n+1} = j , X_n = i)\\ &= \sum_{i \in S} P(X_{n+1} = j \mid X_n = i) P(X_n = i) = \sum_{i \in S} \pi_i^{(n)} p_{ij} \end{align}となるので、ひとまとめにして
となる。これにより、マルコフ連鎖の時間発展は、遷移確率行列を用いて簡単に計算できる。
今後の議論を進めるため、いくつか言葉を定義しよう。まず、状態 $i$ から $j$ へ、また状態 $j$ から $i$ へ有限の時間で遷移できるとき、状態 $i$ と $j$ は相互到達可能であるという。
すべての状態が互いに相互到達可能であるとき、そのようなマルコフ連鎖は既約であるという。
次に、状態 $i$ から出発して $i$ 自身へ戻る時刻 $k$ を考えよう。戻れる時刻 $k$ の最大公約数を、状態 $i$ の周期 $d(i)$ と定義する。
ただし、上の周期は戻れる時刻が存在する場合に定義する。$d(i)>1$ なら、状態 $i$ は周期的である。$d(i)=1$ なら、状態 $i$ は非周期的である。例えば $p_{ii}>0$ なら、状態 $i$ は（ $1$ 回で戻れるので）非周期的である。
このときまず、以下の命題が成り立つ。
命題： 状態 $i,j$ が相互到達可能ならば、 $i,j$ は同じ周期をもつ（つまり $d(i)=d(j)$ ）。
証明しよう。まず仮定より、ある時刻 $n,m \geq 1$ が存在して、$P(X_n=j \mid X_0=i)>0$ および $P(X_m=i \mid X_0=j)>0$ が成り立つ。
状態 $i$ から $j$ を経由して再び $i$ に戻る経路を考えると、$n+m$ は状態 $i$ の再帰時刻であるから、$d(i)$ は $n+m$ の約数である。
さらに、状態 $j$ から $j$ に戻る任意の再帰時刻を $k>0$ とする。このとき、状態 $i$ から $j$ に移動し、$j$ に戻る経路を経て、再び $i$ に戻ることができるため、$n+k+m$ も状態 $i$ の再帰時刻となる。したがって、$d(i)$ は $n+k+m$ の約数でもある。
以上から差を取ることにより、$d(i)$ は $k$ の約数である。これは状態 $j$ の任意の再帰時刻 $k$ について成立するため、$d(i)$ は $d(j)$ の約数となる。同様に、$d(j)$ は $d(i)$ の約数となることから、$d(i)=d(j)$ が成り立つ。
以上から、次の命題が成り立つ。
命題： 既約なマルコフ連鎖では、すべての状態の周期が等しい。
次に、状態 $i$ から出発したとき、$n$ 回目に初めて状態 $j$ に到達する確率を $f_{ij}^{(n)}$ と書く。すなわち、
\begin{align} f_{ij}^{(n)} &= P( X_n =j,\ X_k \neq j ,\ k=1,2,\dots,n-1 \mid X_0=i) \end{align}である。これに基づいて、
\begin{align} f_{ij} &= \sum_{n=1}^{\infty} f_{ij}^{(n)} \end{align}とかく。これは状態 $i$ を出発して、いつか状態 $j$ へ到達する確率を表す。したがって、$f_{ii}$ は状態 $i$ から出発していつか戻ってくる確率を表す。$f_{ii} = 1$ なら状態 $i$ は再帰的であるといい、$f_{ii} < 1$ なら状態 $i$ は一時的であるという。
次に、$i$ から $j$ への平均到達時間 $\mu_{ij}$ を考える。初めて状態 $j$ に到達する正の時刻を $T_j$ とし、到達しない場合は $T_j=\infty$ とおく。$\mu_{ij}=E[T_j\mid X_0=i]$ と定義すると、到達確率 $f_{ij}=1$ の場合には、
が成り立つ。一方、$f_{ij}<1$ の場合は、到達しない確率が正であるため $\mu_{ij}=\infty$ となり、上の和だけでは期待値を表せない。特に $i=j$ の場合、$\mu_{jj}$ は平均帰還時間である。状態 $j$ が再帰的であり、かつ $\mu_{jj} < \infty$ なら正再帰的、$\mu_{jj}=\infty$ なら零再帰的であるといわれる。
正再帰性については、以下の命題が成り立つ。
命題： 状態空間 $S$ が有限集合・既約なマルコフ連鎖では、すべての状態 $i \in S$ が正再帰的である。
ようは、有限個の状態がすべてつながっているなら、各状態 $i$ には確率 $1$ で戻り、その平均帰還時間 $\mu_{ii}$ も有限になる。この結論に非周期性は必要ない。たとえば、$2$ 状態を交互に行き来する連鎖は周期 $2$ だが、どちらの状態にも必ず $2$ ステップで戻るので正再帰的である。
では、十分長く動かした後に各状態にいる確率はどうなるだろうか。以下、 $p_{ij}^{(n)}=P(X_n=j\mid X_0=i)$ と書く。
命題： 既約・再帰的・非周期的なマルコフ連鎖では、任意の状態 $i,j$ について次式が成り立つ。ただし $\mu_{jj}=\infty$ のときは $1/\mu_{jj}=0$ とする。
\begin{align} \lim_{n \to \infty} p_{ij}^{(n)} = \frac{1}{\mu_{jj}} \label{eq:lim_pij} \end{align}\eqref{eq:lim_pij} の右辺が到着点 $j$ のみに依ることに注目しよう。十分長い時間では「どこから出発したか」の影響が薄れ、状態 $j$ のみによって値が決まる。平均帰還時間 $\mu_{jj}$ が短い状態ほど、長い目で見た滞在割合は大きくなる。非周期性は、時刻の偶奇などに応じた振動が残ることを防ぐために必要な条件である。
ただし、零再帰的な場合には、どの固定状態 $j$ についても極限は $0$ になる。これは確率が消えるという意味ではなく、可算無限の状態空間へ広く分散し、各状態の確率が小さくなるということである。この場合、確率の総和が $1$ の定常分布は存在しない。
そこで、次に既約・正再帰的な連鎖を考えよう。この場合には、平均帰還時間の逆数が定常分布の各成分を与える。
命題： 既約・正再帰的なマルコフ連鎖には、一意な定常分布 $\boldsymbol\pi$ が存在し、任意の状態 $j\in S$ について次式が成り立つ。非周期性は必要ない。
\begin{align} \boldsymbol\pi P&=\boldsymbol\pi, \quad \sum_{j\in S}\pi_j=1, \quad \pi_j=\frac{1}{\mu_{jj}}>0 \label{eq:pi_j} \end{align}定常分布から出発すると、個々の軌跡では状態が変わったり、同じ状態に留まったりするが、全体としての確率分布は $1$ ステップ進んでも変化しない。また、平均帰還時間が短い状態ほど、定常確率が大きい。
なお、実際に $\mu_{jj}$ を計算したいときは、 $\mu_{jj}$ を直接求めるより上の \eqref{eq:pi_j} を連立して解いて、$\pi_j$ の逆数から求める方が簡単である点にも注意しておこう。
命題： 既約・正再帰的・非周期的なマルコフ連鎖では、任意の状態 $i,j$ について次式が成り立つ。特に有限状態の場合は、既約性と非周期性を仮定すればよい。
\begin{align} \lim_{n \to \infty} p_{ij}^{(n)} = \pi_j \label{eq:lim_pij2} \end{align}\eqref{eq:lim_pij2} は、どの状態から出発しても、十分時間がたつと同じ定常分布に近づくという結論である。
一方で、周期的なマルコフ連鎖では上の命題は成り立たないが、時間平均を取ることで定常分布への収束が保証される。
命題： 既約・正再帰的なマルコフ連鎖では、任意の状態 $i,j$ について次式が成り立つ。非周期性は必要なく、特に有限状態の場合は、既約性だけを仮定すればよい。
\begin{align} \lim_{n \to \infty} \frac{1}{n+1} \sum_{k=0}^{n} p_{ij}^{(k)} = \pi_j \label{eq:lim_pij3} \end{align}最後に、マルコフ連鎖の可逆性を定義しておこう。確率分布 $\boldsymbol\pi=(\pi_i)_{i\in S}$ と遷移確率行列 $P=(p_{ij})$ が、すべての状態 $i,j$ について次の詳細釣り合い条件を満たすとき、$P$ は $\boldsymbol\pi$ に関して可逆であるという。
なお、上の式では添字の重複があるが、和を取っていないことに注意しよう（アインシュタインの縮約は用いていない）。
実際、詳細釣り合い条件の両辺を $i$ について足し合わせると、
\begin{align} \sum_{i\in S}\pi_i p_{ij} =\pi_j\sum_{i\in S}p_{ji} =\pi_j \end{align}となる。最後の等号では、状態 $j$ からの遷移確率の和が $1$ であることを使った。左辺は $1$ ステップ後に状態 $j$ にいる確率なので、分布が変わらないことが分かる。
命題： $P$ が確率分布 $\boldsymbol\pi$ に関して可逆なら、$\boldsymbol\pi$ は $P$ の定常分布である。
\begin{align} \boldsymbol\pi P=\boldsymbol\pi \end{align}この事実を使うと、定常分布の候補について、状態の組ごとの詳細釣り合い条件を確認すれば定常性を示せる。また、望む分布 $\boldsymbol\pi$ に対してこの条件を満たすように $P$ を設計することで、その分布を定常分布にもつ連鎖を作れる。これは MCMC で用いられる基本的な考え方である。
その他、確率過程にはランダム・ウォークなどの重要な事項があるが、これ以上は専門の教科書に説明を譲る。