em算法 三硬币模型

三硬币模型

设有3枚硬币记作a,b,c。出现正面的概率分别为π,p,q。先掷硬币a，若a为正面选硬币b，反面选c。

E步

观测数据y来自硬币b的概率

M步

$$p=\frac{m{\mu }_1}{m{\mu }_1+(n-m){\mu }_2}$$
$$q=\frac{m(1-{\mu }_1)}{m(1-{\mu }_1)+(n-m)(1-{\mu }_2)}$$
此时计算下一步的${\mu }_1^{(1)}=\frac{\pi p}{\pi p+(1-\pi )q}$
带入$\pi$,p,q后：
$${\mu }_1^{(1)}=\frac{\frac{m}{n}{\mu }_1}{\frac{m}{n}{\mu }_1+(1-\frac{1}{n}[m{\mu }_1+(n-m){\mu }_2])\frac{m(1-{\mu }_1)}{m(1-{\mu }_1)+(n-m)(1-{\mu }_2)}}$$
$${\mu }_1^{(1)}=\frac{{\mu }_1}{{\mu }_1+(n-[m{\mu }_1+(n-m){\mu }_2])\frac{1-{\mu }_1}{m(1-{\mu }_1)+(n-m)(1-{\mu }_2)}}$$
$${\mu }_1^{(1)}=\frac{{\mu }_1}{{\mu }_1+(n-n{\mu }_2+m{\mu }_2-m{\mu }_1)\frac{1-{\mu }_1}{m-m+n-n{\mu }_2+m{\mu }_2-m{\mu }_1}}$$
$${\mu }_1^{(1)}=\frac{{\mu }_1}{{\mu }_1+1-{\mu }_1}$$
$${\mu }_1^{(1)}={\mu }_1$$

个人觉得书上三硬币的例子并不是迭代算法，第二次迭代生成的${\mu }_1$等于第一次生成的${\mu }_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;
} 

使用代码模拟十次迭代