EM算法例子 matlab

EM算法

EM算法 实例讲解

%{
实现的是一个AB硬币的例子
thetaA和thetaB表示硬币A和B正面向上的概率
A的计算结果是在每一轮投掷时,选择的是硬币A的概率
%}
clear
clc
data = [[5 5];[9 1];[8 2];[4 6];[7 3]];
m_step = 10;
PA = zeros(1,5);
PB = zeros(1,5);
A = zeros(5,1);
thetaA_start = 0.5;
thetaB_start = 0.4;
resultA = zeros(1,m_step);
resultB = zeros(1,m_step);
thetaA = thetaA_start;
thetaB = thetaB_start;

for iteration = 1:m_step
    for k = 1:length(data)
        onedata = data(k,:);
        PA(k) = thetaA^onedata(1)*(1-thetaA)^onedata(2);
        PB(k) = thetaB^onedata(1)*(1-thetaB)^onedata(2);
        A(k) = PA(k)/(PA(k)+PB(k));
    end
    AH = sum(A.*data(:,1));
    AT = sum(A.*data(:,2));
    thetaA = AH/(AH+AT);
    resultA(iteration) = thetaA;
    
    B = 1-A;
    BH = sum(B.*data(:,1));
    BT = sum(B.*data(:,2));
    thetaB = BH/(BH+BT);
    resultB(iteration) = thetaB;
end
plot(1:m_step,resultA)
hold on
plot(1:m_step,resultB)

A
resultA(m_step)
resultB(m_step)

结果图
如何感性地理解EM算法?

  • 4
    点赞
  • 37
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
EM算法(Expectation-Maximization Algorithm)是一种迭代算法,用于在具有潜变量(latent variables)的统计模型中进行参数估计。它常用于无监督学习和模型的最大似然估计。可以用MATLAB编写实现。下面是一个简单的例子: 假设我们观测到一组数据X={x1,x2,...,xn},我们知道这些数据是从一个高斯混合模型中生成的,但是我们不知道这个混合模型的参数,即每个高斯分布的均值和方差,以及每个高斯分布的权重。我们可以使用EM算法来估计这些参数。 步骤如下: 1.初始化模型参数,包括每个高斯分布的均值、方差和权重。 2.在E步骤中,计算每个数据点属于每个高斯分布的概率。这可以通过计算后验概率来完成,即对于每个数据点和每个高斯分布,计算P(该数据点属于该高斯分布|模型参数)。这可以使用贝叶斯公式来计算。 3.在M步骤中,使用上一步计算出的后验概率来更新模型参数。具体来说,更新每个高斯分布的均值、方差和权重,使得对数似然函数最大化。 4.重复步骤2和步骤3,直到模型参数达到收敛。 下面是MATLAB代码示例: % 生成一组数据 mu1 = [1 2]; sigma1 = [1 0.5;0.5 2]; mu2 = [-3 -5]; sigma2 = [2 0;0 1]; X = [mvnrnd(mu1,sigma1,100);mvnrnd(mu2,sigma2,100)]; % 初始化模型参数 k = 2; % 高斯分布的个数 pi = ones(1,k)/k; % 权重 mu = [ones(k,1)*mean(X);ones(k,1)*mean(X)+2]; % 均值 sigma = repmat(eye(2),1,1,k); % 方差 % EM算法迭代 for iter = 1:10 % 设定最大迭代次数 % E步骤 for i = 1:k pdf(:,i) = mvnpdf(X,mu(:,i)',sigma(:,:,i)); % 计算每个数据点属于每个高斯分布的概率 end gamma = repmat(pi,size(X,1),1).*pdf; % 计算后验概率 gamma = gamma./repmat(sum(gamma,2),1,k); % 归一化 % M步骤 Nk = sum(gamma,1); % 每个高斯分布的数据点数 pi = Nk/size(X,1); % 更新权重 for i = 1:k mu(:,i) = sum(repmat(gamma(:,i),1,size(X,2)).*X)/Nk(i); % 更新均值 X_centered = X-repmat(mu(:,i)',size(X,1),1); sigma(:,:,i) = (X_centered.*(repmat(gamma(:,i),1,size(X,2)).*X_centered))'/Nk(i); % 更新方差 end end % 绘制结果 figure; scatter(X(:,1),X(:,2),'filled'); hold on; ezcontour(@(x,y)pdf_2d(x,y,mu,sigma,pi),[-10 10],[-10 10]); title('Gaussian Mixture Model'); function p = pdf_2d(x,y,mu,sigma,pi) p = zeros(size(x)); for i = 1:size(mu,2) p = p+pi(i)*mvnpdf([x y],mu(:,i)',sigma(:,:,i)); end end
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值