2019 人文　（統計検定1級応用　過去問の解説です）

1. 条件付期待値・分散
2. クラスター分析
3. 項目反応理論
4. 分散構造分析
5. カイ二乗適合度検定
   　

---

2019人文① 条件付期待値・分散

２変量正規分布と条件付期待値の理解が肝、理論的背景の理解を求められる。

［条件付期待値］ 　（大学生の確率統計 p130参照)

Xの確率密度関数が、 $f_{X}( x)$ であるとき、

　　　　　$a\leq X\leq b$ における条件付き期待値は、 　　　$E\left( X| a\leq X\leq b\right) =\displaystyle \int _{a}^{b}x\cdot \dfrac{f_{X}\left( x\right) dx}{P\left( a\leq X\leq b\right) }$

　　　　　条件付き確率密度関数は、 　　　　　　　　　　　　　　　$f_{X| a\leq x\leq b}\left( x\right) =\dfrac{f_{X}\left( x\right) }{P\left( a\leq X\leq b\right) }$

　積分を $P\left( a\leq X\leq b\right)$ で割っているだけなので、期待値や分散の線型性は成立する。

［参考］
連続型の $E\left( Y| X\right)$ のとき、 $\displaystyle E\left(Y| X=x\right) =\int ^{\infty }_{-\infty }yf_{Y|X}\left(y| x\right) =\int ^{\infty }_{-\infty }y\dfrac{f_{(X,Y)}\left( x,y\right) }{f_{X}\left( x\right) }dy$　　　これは後で説明する。

数学ノート 様 – ２変量正規分布の条件付き期待値,分散
　　　　　 https://math-note.com/conditional-probability-of-multivariate-normal-distribution/

　

(1) 受験生の得点 $X\sim N\left( 100,20^{2}\right)$ であるので、$Z=\dfrac{X-100}{20}\sim N(0,1)$

Zで計算し、 $X=20Z+100$ によってXに戻せばよい。

偏差値は、正規分布N(50,10)であるから、

　　　　　 A君のZ値は、 $Z=0.4$ 　　　　　得点は $X=20\times 0.4+100=108$

A君の位置は、全体では上位　 $P(Z\geq 0.4)=0.3446$ に位置する。

合格者の割合は、 $P(Z\geq 0)=0.5$ だから、合格者の中では、 $\frac{0.3446}{0.5}=0.6892$ に位置する。

入学者の割合は、0.4 だから、入学者の中では、 $\dfrac{0.3446-1}{0.4}=0.6115$ に位置する。

(2) 解答参照

(3) Xの分布関数は、$F(x)=P\left( X\leq x\right) =P\left( Z\leq \dfrac{x-100}{20}\right) =\displaystyle \int ^{\dfrac{x-100}{20}}_{-\infty }\varphi\left( z\right) dz$

確率密度関数は、$f_X(x)=F'(x)=\frac{1}{20}\varphi\left(\dfrac{x-100}{20}\right)$

$X\geq 100$ での確率密度関数は、

　　　　　$f_{X| X\geq 100}\left( x\right) =\dfrac{f_{X}\left( x\right) }{P\left(X\geq 100\right) } =\frac{1}{20}\varphi\left(\dfrac{x-100}{20}\right) \div \frac{1}{2} =\frac{1}{10}\varphi\left(\dfrac{x-100}{20}\right)$

(4) 条件付き期待値は、

　　　　　$\begin{split} E(X|X\geq 100)&=\int ^{\infty }_{100}x\dfrac{1}{10}\varphi\left( \dfrac{x-100}{20}\right) dx\\ &=\int ^{\infty }_{0}(20z+100)\dfrac{1}{10}\varphi\left(z\right) 20dz\\ &=2\int ^{\infty }_{0}(20z+100)\varphi\left(z\right)dz =40\int ^{\infty }_{0}z\varphi\left(z\right)dz+100 \end{split}$

で計算してもよいが、Zを利用するのがよい。即ち、

　　　　　$E(X|X\geq 100)=E(20Z+100|Z\geq 0) =20\cdot E(Z|Z\geq 0)+100$

　　　　　$E(Z|Z\geq0) =\int ^{\infty }_{0}z\frac{\varphi\left(z\right)}{P(Z\geq 0)}dz =2\int ^{\infty }_{0}z\varphi\left(z\right)dz$

ここで、 $\left(exp\left({-\dfrac{z^{2}}{2}}\right)\right)'=-z\cdot exp \left({-\dfrac{z^{2}}{2}}\right)$ より、$\displaystyle \int z\varphi\left( z\right) =-\varphi \left( z\right)$

　　　　　$E(Z|Z\geq0)=2\left[-\varphi \left( z\right) \right]_{0}^{\infty} =2\varphi(0)=2\frac{1}{\sqrt{2\pi}}=\sqrt{\frac{2}{\pi}}$

よって、　$E(X|X\geq 100)=20\sqrt{\frac{2}{\pi}}+100=115.96$

条件付き分散は、 $X=20Z+100$ より、

　　　　　$V(X|X\geq100)=V(20Z|Z\geq0)=20^2V(Z|Z\geq0)=20^2\left[E(Z^2|Z\geq0)-E(Z|Z\geq0)^2\right]$

　　　　　$\begin{split} E(Z^2|Z\geq0) &=\int ^{\infty }_{0}z^2\frac{\varphi\left(z\right)}{P(Z\geq 0)}dz =2\int ^{\infty }_{0}z^2\varphi\left(z\right)dz\\ &=2\int ^{\infty }_{0}z\cdot z\varphi\left(z\right)dz =2[z(-\varphi(z)]_{0}^{\infty}+2\int_{0}^{\infty}\varphi(z)dz=1 \end{split}$

よって、　$V(X|X\geq100)=20^2\left(1-\frac{2}{\pi}\right)$

［注意］

$E(X^2|X\geq100)=\int ^{\infty }_{0}(20z+100)^2\dfrac{1}{10}\varphi\left(z\right) 20dz =2\int ^{\infty }_{0}(20z+100)^2\varphi\left(z\right)dz =400+4000\sqrt{\frac{2}{\pi}}+10000$

$E(X|X\geq 100)^2=\left(20\sqrt{\frac{2}{\pi}}+100\right)^2=400 \frac{2}{\pi}+4000\sqrt{\frac{2}{\pi}}+10000$

$V(X|X\geq100)=E(X^2|X\geq100)-E(X|X\geq 100)^2=20^2\left(1-\frac{2}{\pi}\right)$　　としてもよいが計算が長い。

条件付き期待値、条件付き分散の性質については、大学生の確率統計p128 、p134、ノート③を参照。

　

2019人文② クラスタリング

Macromill 様 – 階層クラスター分析
　　　　　https://www.macromill.com/service/data_analysis/hierarchical-clustering.html

　　　　　最短距離法を用いてデンドログラムを手書きするとよい。

同志社大学 金明哲様 – Rとクラスター分析
　　　　　 https://www1.doshisha.ac.jp/~mjin/R/Chap_28/28.html　　(Rによる実装 hclust)

（ 最も距離が近い2つの個体間の距離をコーフェン距離とした距離行列がコーフェン行列 ）

　　　　　データ行列　⇒　距離行列　⇒　コーフェン行列　⇒　デンドログラム

(1) 解説参照

Rの出力から、デンドログラム作成の手順を考える。

> man.d<- dist(man);man.d 距離行列

	A	B	C	D	E
B	1.414214
C	2.236068	2.236068
D	8.485281	8.602325	6.403124
E	8.062258	7.81025	5.830952	2.236068
F	10	9.899495	7.810250	2.000000	2.236068

> man.hc$merge クラスター形成の履歴　　　　　

	[,1]	[,2]
[1,]	-1	-2	①=(AB)	ABを１番目のクラスターに
[2,]	-4	-6	②=(DF)	DFを２番目のクラスターに
[3,]	-3	1	③=(C,①)	Cと①=(AB)を３番目のクラスターに
[4,]	-5	2	④=(E,②)	Eと②=(DF)を４番目のクラスターに
[5,]	3	4	⑤=(③,④)	③と④を５番目のクラスターに

階層的クラスターを形成する階層の情報である。マイナス符号が付いているのが個体の番号、

マイナスが付いていないのがクラスターの番号である。行番号がクラスター形成の順番である。

> man.hc$height クラスター結合の際の距離

[1]

1.414214

2.000000

2.236068

2.236068

5.830952

> man.hc$order デンドログラムの書き順

[1] 3 1 2 5 4 6 　CABEDFの順に書く

> cophenetic(man.hc) #コーフェン行列

	A	B	C	D	E
B	1.414214
C	2.236068	2.236068
D	5.830952	5.830952	5.830952
E	5.830952	5.830952	5.830952	2.236068
F	5.830952	5.830952	5.830952	2.000000	2.236068

コーフェン行列は、属するクラスター間の距離を表す。

実用としては、Rのhclustを利用するのがよい。

［参考］

コーフェン相関係数とは、解析前の距離行列とコーフェン行列の相関を求めたもの

相関係数が高ければ、距離行列と用いた方法のコーフェン行列との歪みが小さい

しかし、必ずしもコーフエン相関係数が高い方法の分類結果が良いとは言えない。

　

(3) 非階層的クラスター分析（kmeans法）

X=(築年数 , 所要時間) をデータとし、ユークリッド距離を用いて分類する。

k-means法なので、２つのクラスターに分けるには、初期代表点を２つ決める必要がある。

ここでは、初期代表点を D(7,9) , E(8,7) と指定されている。

D、E からのユークリッド距離なので、境界線はDEの垂直二等分線である。これにより、

１回目の更新で、DF 、ABCE と分かれる。

各クラスターの座標の平均が代表点なので、

　　　　　DFの代表点は　　　　　$\left(\frac{7+9}{2}\ ,\ \frac{9+9}{2}\right)=(8,9)$

　　　　　ABCEの代表点は 　　 $\left(\frac{1+2+3+8}{4}\ ,\ \frac{3+2+4+7}{4}\right)=(3.5 ,\ 4)$

2回目の更新は、(8,9) (3.5,4) の垂直二等分線

　　　　　$y=-0.9x+7.15$　　によって、ABC、DEFと分かれる。

［参考］

Albert 様　『クラスター分析の手法③（非階層クラスター分析）』を参照
　　　https://www.albert2005.co.jp/knowledge/data_mining/cluster/non-hierarchical_clustering

(4) の k-means法の初期値依存についても、上記サイトに解説がある。

k-means法の一つの短所として、初期値（初期に選択される「核」となるk個のサンプル）依存性があります。同じデータを距離などを同じ条件にして計算しても、初期値が異なるだけで、結果が大きく違うことがあります。従って、よいクラスターを得るためには、初期値を変えて何回か分析を実施し、平均クラスター内距離が最小になる初期値を選択するなど、最適初期値での結果を採用することが望ましいといえます。

Rでの実装については、

@Haruka-Ogawa 様　『Rでクラスタリングしてみた』を参照
　　　https://qiita.com/Haruka-Ogawa/items/fcda36cc9060ba851225

　　　irisデータを用いた実装が紹介されている。

［参考］

判別分析：分類のわかっているデータ（教師データ）があり、データの属性値に基づいて、この分類を判別する手法である。機械学習の分野では、教師あり学習（supervised learning）と呼ばれている。フィッシャーの線形判別など、群間変動を最大化する方法については別途。

（ノート③ 2019 の中の、「涌井」、「機械学習に強くなりたいブログ」など参照）

クラスター分析：分類がわかっていないので、教師無し学習（unsupervised learning）と呼ばれる。

線形判別の理論的説明については、

機械学習に強くなりたいブログ 様
　　　https://www.iwanttobeacat.com/entry/2018/05/13/231628 　参照

　

2019人文③ 項目反応理論

(1) 項目反応関数 　　$P_j(\theta) =\dfrac{1}{ 1+\exp \left[ -a_j(\theta -b_{j}) \right]}$　より　 $P_j(0) =\dfrac{1}{ 1+\exp a_jb_{j} }$　　だから、

$a_jb_j$ が等しい項目は、項目２と項目３である。

(2) 困難度 $b_j$ が最も小さいものは、項目４である。

(3) 識別力 $a_j$ は、項目特性曲線の傾きに対応する。

$1-P_j(\theta) =\dfrac{\exp \left[ -a_j(\theta -b_{j}) \right]}{ 1+\exp \left[ -a_j(\theta -b_{j}) \right]}$ に注意して、

　　　　　$P_j'(\theta) =\dfrac{a_j\exp \left[ -a_j(\theta -b_{j})\right]}{ (1+\exp \left[ -a_j(\theta -b_{j}) \right])^2}=a_jP_j(\theta)(1-P_j(\theta))$

$a_j$ が大きくなれば $P_j'(\theta)$ も大きくなりそうである。

実際、グラフの中心 $P_j(\theta)=0.5$ である点において、$P_j'(\theta) =\frac{a_j}{4}$ となる。

(4) 一般的定義を確認

項目情報関数（Item Information Function）

　　　　　$f_i(\theta)=a_i^2P_i(\theta)(1-P_i(\theta))$　　問題の情報量を表す。 $P_i(\theta)=\dfrac{1}{2}$ で最大となる。

（ラッシュモデルの場合には、項目情報関数は単純に正しい応答確率と不正確な応答確率の積）

姫路獨協大学医療保健学部理学療法学科 様
　　　　　　　　　　 http://kinki59.sakura.ne.jp/download/P-101_C000131.pdf　　　によれば、

　　　　　$f_j(\theta)=\dfrac{P_j'(\theta) ^2}{P_j(\theta)(1-P_j(\theta))}$　を定義としているが、これは

　$P_j'(\theta)^2 =a_j^2P_j(\theta)^2(1-P_j(\theta))^2$　より、 $\dfrac{P_j'(\theta) ^2}{P_j(\theta)(1-P_j(\theta))}=a_j^2P_j(\theta)(1-P_j(\theta))$ であるから、

右辺を定義としても、左辺を定義としてもよいということである。

（３パラメーターの場合、式は変わる）

$f_j(\theta)$ は問題による $\theta$ の推定精度を表す。受験生の能力値 $\theta$ に応じて値が変わることに注意。

テスト情報関数（Test Information Function）

　　　　　$\displaystyle I(\theta)=\sum_{i=1}^n f_i(\theta)$　　　　　テスト全体の情報量を表す。

（受験生の能力値のレベルに応じて値が変わることが大きな特徴である）

(5) 　$c_j$ は当て推量パラメータと呼ばれる。（解答参照）

　

2019人文④ 共分散構造分析

標準解とは、観測変数および潜在変数を標準化した解であり、因果の強さを見るのに適する。

(1) パス図は変数同士の関係を表している。

楕円で囲まれた変数は潜在変数、四角に囲まれた変数は観測変数、円で囲まれた変数は誤差を表す。パス図からそれらの関係を読み取る。

　　　　　観測変数は、 $y_1=0.652\ f_1+e_1$　　　　誤差分散　 $V(e_1)=0.572$

潜在変数は標準化されているから、

　　　　　$V(y_1)=(0.654)^2+0.572=1$

観測変数が標準化されているパス図を標準解という。

　　

(2) 潜在変数どうしの関係は、

　　　　　$f_2=0.640\ f_1+e_{10}$　　　 $f_3=0.014\ f_1+0.032\ f_2+e_{11}$

　

(3)(4) 構造方程式　 $f=Af+e$ は、

　　　　　$\begin{pmatrix} f_1 \\ f_2 \\ f_3 \end{pmatrix} =\begin{pmatrix} 0 & 0 & 0\\ 0.640 & 0 & 0\\ 0.014 & 0.032 & 0 \end{pmatrix} \begin{pmatrix} f_1 \\ f_2 \\ f_3 \end{pmatrix} +\begin{pmatrix} f_1\\ e_{10} \\ e_{11} \end{pmatrix}$　　　（ 　$f_1$ はランダムに変化するので誤差に含める ）

(5) 同値モデルとは、モデルが異なっても分散共分散行列が一致するモデルである。

このような場合は適合度が同じになる。（解答参照）

　　　

2019人文⑤ カイ二乗適合度検定

AVILEN AI Trend 様– 適合度検定をカイ二乗検定で行う例題と手順
　　　　　https://ai-trend.jp/basic-study/hypothesis-testing/goodness-of-fit-test-by-chi-test/

統計科学研究所 様 – 適合度検定の理論的説明
　　　　　https://statistics.co.jp/reference/StatisticalDataScienceA/Chapter7nonpassword.pdf　　　　　

適合度検定とは観測度数が理論比率に基づいて得られているかを検証する仮説検定である。

カテゴリナンバー $i$ 番目の観測度数を $n_i$ 期待度数を $E_i$ とするとき、

X	1	2	‥	r
観測度数	$n_1$	$n_2$	‥	$n_r$
期待度数	$E_1$	$E_2$	‥	$E_r$

カイ二乗値 　　$\chi^2 = \sum^{r}_{i=1}\frac{(n_{i}-E_{i})^2}{E_{i}}\sim \chi^2(r-1)$　　　 自由度 r-1 のカイ二乗分布に従う

(1) 　カイ二乗値は、4.0 となる。（解答参照）

(2) 　度数がc倍になると、カイ二乗値もc倍になる。

　　　　　$\chi^2 = \sum^{r}_{i=1}\frac{(cn_{i}-cE_{i})^2}{cE_{i}}=c \sum^{r}_{i=1}\frac{(n_{i}-E_{i})^2}{E_{i}}$

$k=1$ のときのカイ二乗値は1 だから、 $k$ のときのカイ二乗値は $k$

$k>7.81$ より、 $K=8$ となる。

(3) 理論的根拠は大阪大学 鈴木穣 教授による説明がある。（難）
　　　　　http://joe.bayesnet.org/wp-content/uploads/2017/05/2017-5-25-1.pdf　

テキストの解答も鈴木先生の方法に沿っている。以下の説明は鈴木先生の解説概略である。

$X=0,1,2,\cdots,(\alpha-1)$ は次の表のような多項分布に従う。（ $p_0+p_1+p_2+\cdots+p_{\alpha-1}=1$ ）

X	0	1	$\cdots$	$j$	$\cdots$	$\alpha-1$
P	$p_0$	$p_1$	$\cdots$	$p_j$	$\cdots$	$p_{\alpha-1}$

n回の試行による $X_1,X_2,\cdots,X_n$ は独立に上記の多項分布に従うとする。

結果 $j$ が起こる回数 $w_j$ を次のように数式化する。

$k$ 回目の試行で結果 $j$ が起きれば１起きなければ０をとる確率変数を

　　　　　$I(X_k=j) =\begin{cases}1\ (X_k=j)\\ 0\ (X_k\neq j)\\ \end{cases}$　　　　とすると、 $I(X_k=j)\sim Ber(p_j)$ 　ベルヌーイ分布に従い、

n 回の試行で $j$ が起きる回数は、

　　　　　$w_j=\sum_{k=1}^{n}I(X_k=j)$　 で表され、 二項分布 $w_j\sim B(n,p_j)$ に従う。

nが大きければ、 $B(n,p_j)\sim N\left(np_j,np_j(1-p_j)\right)$ 正規近似できる。即ち、

　　　　　$\dfrac{w_j-np_j}{\sqrt{np_j(1-p_j)}}\rightarrow N(0,1)$　　　よって、$Z_j=\dfrac{w_j-np_j}{\sqrt{np_j}}\rightarrow N\left(0,(1-p_j)\right)$

X	0	1	$\cdots$	$j$	$\cdots$	$\alpha-1$
実現回数	$w_0$	$w_1$	$\cdots$	$w_j$	$\cdots$	$w_{\alpha-1}$
期待度数	$np_0$	$np_1$	$\cdots$	$np_j$	$\cdots$	$np_{\alpha-1}$

　　　　$\chi^2=\sum_{j=0}^{\alpha-1}\left(\dfrac{w_j-np_j}{\sqrt{np_j}}\right)^2=\sum_{j=0}^{\alpha-1}Z_j^2$　　（α個の正規分布の平方和）

ある直交変換　 $(Z_j)\rightarrow (V_j)$ によって、（ここは難しいので略）

　　　　　$\chi^2=\sum_{j=1}^{\alpha-1}V_j^2$　　（ $V_j$ は独立な標準正規分布）

$(\alpha-1)$ 個の独立な標準正規分布の平方和になるので、自由度 $(\alpha-1)$ のカイ二乗分布になる。

［参考］
AVILen AItrend 様 「多項分布とは？期待値・分散・共分散」
　　　　　　　　　　https://ai-trend.jp/basic-study/multinomial-distribution/m-parameter-derivation/

多項分布の分散共分散行列

　　　　　$E(X_{i})= \displaystyle np_{i}~ ~,~ Var(X_{i}) = \displaystyle np_{i}(1 - p_{i}) ~~ (i = 1, …, k)$　　　　$Cov(X_{i}, X_{j}) =\displaystyle -np_{i}p_{j} ~~ (i < j)$

難しいとした直交変換はヘルマート変換と呼ばれる。（当サイト「正規分布と関連分布」参照）

　

(4) （血液型の遺伝理論については専門の書籍にて確認ください）

血液型遺伝子の存在比率は、 $O\rightarrow r,A\rightarrow p,B\rightarrow r$ とし、これを推定することを目標とする。

実際に観測される人の血液型は、２つの遺伝子配列によって決まり、

　　　　　O型＝OO、A型＝AA＋(AO+OA)、B型=BB＋(BO＋OB)、AB型＝AB＋BA

２つの遺伝子配列の度数を、

　　　　　$O\rightarrow f_{OO}\ ,\ A\rightarrow f_{AA}+f_{AO}\ ,\ B\rightarrow f_{BB}+f_{BO}\ ,\ AB\rightarrow f_{AB}$　　とする。

血液型の観測度数を $n_O,n_A,n_B,n_{AB}$ 、遺伝子型の度数を $f_{OO},f_{AA},f_{AO},f_{BB},f_{BO},f_{AB}$

　　　　　$n_O=f_{OO}=Nr^2$　　　　　$n_A=f_{AA}+f_{AO}=Np^2+2Npr$

　　　　　$n_{AB}=f_{AB}=2Npq$　　　 $n_B=f_{BB}+f_{BO}=Nq^2+2Nqr$

　　　　　（ $f_{AO}$ はAO型とOA型の合計であることに注意 ）

　　　　　$n_O+n_A+n_B+n_{AB}=f_{OO}+(f_{AA}+f_{AO})+(f_{BB}+f_{BO})+f_{AB}=N$

$n_O,n_A,n_B,n_{AB}$ は観測される。

$f_{OO}\ ,\ f_{AB}$ は観測値がそのまま予測値である。

$f_{AA},f_{AO},f_{BB},f_{BO}$ は観測されないので推定が必要である。

　

(ⅰ) 血液型の観測度数 $n_O,n_A,n_B,n_{AB}$ が得られたとき、

　　　　　遺伝子型の度数 $f_{OO},f_{AA},f_{AO},f_{BB},f_{BO},f_{AB}$ から、母数 $r,p,q$ を推定する。

尤度関数を $L(r,p,q)$ とする。（解説参照）

条件 $r+p+q-1=0$ のもとで、対数尤度関数　 $\log L(r,p,q)$ を最大にする $r,p,q$ を求めるには、

ラグランジュの未定乗数法を用いる。

ラグランジュ関数は $Q=\log L(r,p,q)-\lambda(r+p+q-1)$

$Q=2f_{OO} \log r+2f_{AA} \log p+f_{AO} \log 2pr+2f_{BB} \log q+f_{BO}\log 2qr+f_{AB}\log 2pq -\lambda(r+p+q-1)+c$

（ $c$ は定数 ）

$Q$ を $r,p,q$ で偏微分することにより、推定値 $\widehat r,\widehat p,\widehat q$ が求まる。

　　　　　$\hat{r}=\dfrac{2f_{OO}+f_{AO}+f_{BO}}{2N}~~,~~\hat{p}=\dfrac{2f_{AA}+f_{AO}+f_{AB}}{2N}~~,~~ \hat{q}=\dfrac{2f_{BB}+f_{BO}+f_{AB}}{2N}$　　　（解答参照）

これは最大化なのでMステップである。

(ⅱ) 血液型の観測度数 $n_O,n_A,n_B,n_{AB}$ が得られたとき、

　　　　　母数 $r,p,q$ から、 $f_{OO},f_{AA},f_{AO},f_{BB},f_{BO},f_{AB}$ を推定する。

　　　　　$E(f_{AA})=Np^2$　　　　　$E(f_{AO})=n_A-Np^2$

　　　　　$E(f_{BB})=Nq^2$　　　　　$E(f_{BO})=n_B-Nq^2$

　　　　　$E(f_{OO})=Nr^2$　　　　　$E(f_{AB})=N\cdot2pq$

これは推定なのでEステップである。

［EMアルゴリズム］

血液型の観測度数 $n_O,n_A,n_B,n_{AB}$ が得られたとき、

$f_{OO},f_{AA},f_{AO},f_{BB},f_{BO},f_{AB}$ と $r,p,q$ の関係を用いてEMアルゴリズムで母数 を推定する。

手順を具体的に書いてみよう。

$f_{OO},f_{AA},f_{AO},f_{BB},f_{BO},f_{AB}$ のうち、 $f_{OO}=n_O\ ,\ f_{AB}=n_{AB}$ とするのが妥当である。

$f_{AA},f_{BB}$ の適当な初期値を定め、 $f_{AA}^{(0)} \ ,\ f_{BB}^{(0)}$ とし、$f_{AO}^{(0)}=n_A-f_{AA}^{(0)}$ 、$f_{BO}^{(0)}=n_B-f_{BB}^{(0)}$ とすると、

(ⅰ)の方法で、

　　　　　遺伝子型の度数 $f_{OO},f_{AA}^{(0)},f_{AO}^{(0)},f_{BB}^{(0)},f_{BO}^{(0)},f_{AB}$ から、母数 $r,p,q$ を推定できる。

母数 $r,p,q$ の推定値は、

　　　　　$r^{(1)}=\dfrac{2f_{OO}+f_{AO}^{(0)}+f_{BO}^{(0)}}{2N}$　 、 $p^{(1)}=\dfrac{2f_{AA}^{(0)}+f_{AO}^{(0)}+f_{AB}}{2N}$　 、$q^{(1)}=\dfrac{2f_{BB}^{(0)}+f_{BO}^{(0)}+f_{AB}}{2N}$

(ⅱ) の方法で、この推定値 $r^{(1)},p^{(1)},q^{(1)}$ から、

　　　　　$f_{AA}^{(1)}=N\cdot \left(p^{(1)}\right)^2$　　　　　$f_{AO}^{(1)}=n_A-f_{AA}^{(1)}=n_A-N\cdot\left(p^{(1)}\right)^2$

　　　　　$f_{BB}^{(1)}=N\cdot \left(q^{(1)}\right)^2$　　　　　$f_{BO}^{(1)}=n_B-f_{BB}^{(1)}=n_B-N\cdot \left(q^{(1)}\right)^2$ 　　

　　　　　$f_{OO}=n_O\ ,\ f_{AB}=n_{AB}$　としてそのまま使用することにする。

同じ操作を繰り返す。

(ⅰ)の方法で、

　　　　　遺伝子型の度数 $f_{OO},f_{AA}^{(1)},f_{AO}^{(1)},f_{BB}^{(1)},f_{BO}^{(1)},f_{AB}$ から、母数 $r,p,q$ を推定できる。

母数 $r,p,q$ の推定値は、

　　　　　$r^{(2)}=\dfrac{2f_{OO}+f_{AO}^{(1)}+f_{BO}^{(1)}}{2N}$ 　、 $p^{(2)}=\dfrac{2f_{AA}^{(1)}+f_{AO}^{(1)}+f_{AB}}{2N}$ 　、 $q^{(2)}=\dfrac{2f_{BB}^{(1)}+f_{BO}^{(1)}+f_{AB}}{2N}$

これを繰り返し、収束するまで進むことにより、母数 $r,p,q$ の推定値を得る。