em算法三硬币模型

最新推荐文章于 2024-04-21 11:36:12 发布

Peiris_

最新推荐文章于 2024-04-21 11:36:12 发布

阅读量781

点赞数 2

文章标签：算法

本文链接：https://blog.csdn.net/Peiris_/article/details/126632071

版权

三硬币模型

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

E步

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

其中y的取值为0或1。
当y取1时 $\mu_1=\dfrac{\pi p}{\pi p +(1-\pi)q}$
当y取0时 $\mu_0=\dfrac{\pi (1-p)}{\pi (1-p) +(1-\pi)(1-q)}$

M步

在这里插入图片描述

假设在n个数据中有m个1，则有n-m个0；
$\pi =\frac1 n [m\mu_1+(n-m)\mu_2]$

$p=\dfrac{m\mu_1}{m\mu_1+(n-m)\mu_2}$
$q=\dfrac{m(1-\mu_1)}{m(1-\mu_1)+(n-m)(1-\mu_2)}$
此时计算下一步的 $\mu_1^{(1)}=\dfrac{\pi p}{\pi p +(1-\pi)q}$
带入 $\pi$ ,p,q后：
$\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)}}$
$\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)}}$
$\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 }}$
$\mu_1^{(1)}=\cfrac{ \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;
}