三硬币模型
设有3枚硬币记作a,b,c。出现正面的概率分别为π,p,q。先掷硬币a,若a为正面选硬币b,反面选c。
E步
观测数据y来自硬币b的概率
其中y的取值为0或1。
当y取1时
μ
1
=
π
p
π
p
+
(
1
−
π
)
q
\mu_1=\dfrac{\pi p}{\pi p +(1-\pi)q}
μ1=πp+(1−π)qπp
当y取0时
μ
0
=
π
(
1
−
p
)
π
(
1
−
p
)
+
(
1
−
π
)
(
1
−
q
)
\mu_0=\dfrac{\pi (1-p)}{\pi (1-p) +(1-\pi)(1-q)}
μ0=π(1−p)+(1−π)(1−q)π(1−p)
M步
假设在n个数据中有m个1,则有n-m个0;
π
=
1
n
[
m
μ
1
+
(
n
−
m
)
μ
2
]
\pi =\frac1 n [m\mu_1+(n-m)\mu_2]
π=n1[mμ1+(n−m)μ2]
p
=
m
μ
1
m
μ
1
+
(
n
−
m
)
μ
2
p=\dfrac{m\mu_1}{m\mu_1+(n-m)\mu_2}
p=mμ1+(n−m)μ2mμ1
q
=
m
(
1
−
μ
1
)
m
(
1
−
μ
1
)
+
(
n
−
m
)
(
1
−
μ
2
)
q=\dfrac{m(1-\mu_1)}{m(1-\mu_1)+(n-m)(1-\mu_2)}
q=m(1−μ1)+(n−m)(1−μ2)m(1−μ1)
此时计算下一步的
μ
1
(
1
)
=
π
p
π
p
+
(
1
−
π
)
q
\mu_1^{(1)}=\dfrac{\pi p}{\pi p +(1-\pi)q}
μ1(1)=πp+(1−π)qπp
带入
π
\pi
π,p,q后:
μ
1
(
1
)
=
m
n
μ
1
m
n
μ
1
+
(
1
−
1
n
[
m
μ
1
+
(
n
−
m
)
μ
2
]
)
m
(
1
−
μ
1
)
m
(
1
−
μ
1
)
+
(
n
−
m
)
(
1
−
μ
2
)
\mu_1^{(1)}=\cfrac{\cfrac m n \mu_1}{\cfrac m n \mu_1+\big( 1-\frac 1 n[m\mu_1+(n-m)\mu_2] \big)\cfrac{m(1-\mu_1)}{m(1-\mu_1)+(n-m)(1-\mu_2)}}
μ1(1)=nmμ1+(1−n1[mμ1+(n−m)μ2])m(1−μ1)+(n−m)(1−μ2)m(1−μ1)nmμ1
μ
1
(
1
)
=
μ
1
μ
1
+
(
n
−
[
m
μ
1
+
(
n
−
m
)
μ
2
]
)
1
−
μ
1
m
(
1
−
μ
1
)
+
(
n
−
m
)
(
1
−
μ
2
)
\mu_1^{(1)}=\cfrac{ \mu_1}{ \mu_1+\big( n-[m\mu_1+(n-m)\mu_2] \big)\cfrac{1-\mu_1}{m(1-\mu_1)+(n-m)(1-\mu_2)}}
μ1(1)=μ1+(n−[mμ1+(n−m)μ2])m(1−μ1)+(n−m)(1−μ2)1−μ1μ1
μ
1
(
1
)
=
μ
1
μ
1
+
(
n
−
n
μ
2
+
m
μ
2
−
m
μ
1
)
1
−
μ
1
m
−
m
+
n
−
n
μ
2
+
m
μ
2
−
m
μ
1
\mu_1^{(1)}=\cfrac{ \mu_1}{ \mu_1+\big(n-n\mu_2+m\mu_2-m\mu_1 \big)\cfrac{1-\mu_1}{m-m+n-n\mu_2+m\mu_2-m\mu_1 }}
μ1(1)=μ1+(n−nμ2+mμ2−mμ1)m−m+n−nμ2+mμ2−mμ11−μ1μ1
μ
1
(
1
)
=
μ
1
μ
1
+
1
−
μ
1
\mu_1^{(1)}=\cfrac{ \mu_1}{\mu_1+1-\mu_1}
μ1(1)=μ1+1−μ1μ1
μ
1
(
1
)
=
μ
1
\mu_1^{(1)}=\mu_1
μ1(1)=μ1
个人觉得书上三硬币的例子并不是迭代算法,第二次迭代生成的 μ 1 \mu_1 μ1等于第一次生成的 μ 1 \mu_1 μ1。所以后续迭代产生的 π \pi πpq都等于第一次计算得到的 π \pi πpq。
#include<bits/stdc++.h>
using namespace std;
const int n=10000;
int l[n]={0};
long double u[n]={0};
int main(){
int seed = time(0);
srand(seed);
for(int i=0;i<n;++i){
if(rand()%10000<4123){
l[i]=rand()%10000<5201;
}else{
l[i]=rand()%10000<7321;
}
}
long double a=0.1231,b=0.2319,c=0.123;
for(int i=0;i<10;++i){
long double sumu=0,sumuy=0,sumnuy=0,sumnu=0;
long double i1=a*b/(a*b+(1-a)*c),i2=a*(1-b)/(a*(1-b)+(1-a)*(1-c));
cout<<"------"<<i1<<' '<<i2<<endl;
for(int j=0;j<n;++j){
if(l[j]){
u[j]=i1;
sumuy+=u[j];
sumnuy+=1-u[j];
}
else{
u[j]=i2;
}
sumu+=u[j];
sumnu+=1-u[j];
}
cout<<"***"<<endl;
cout<<sumu<<endl;
cout<<sumuy<<endl;
cout<<sumnuy<<endl;
cout<<sumnu<<endl;
cout<<"***"<<endl;
a=sumu/n;
b=sumuy/sumu;
c=sumnuy/sumnu;
cout<<a<<"\t"<<b<<"\t"<<c<<"\t"<<endl;
}
return 0;
}
使用代码模拟十次迭代
可以看到均与第一轮保持一致