EM算法matlab，或者是算法的程序（只要输入data数据便能出结果的哪种，）或者是写出具体的运行步骤。3q

x=VarName2'

length=size(x,2)

pi0=0.4

p0=0.6

q0=0.7

pi=zeros(1,length)

a=zeros(1,length) %mu和xi的乘积

b=zeros(1,length) %1-mu和xi的乘积

p1=zeros(1,length)

p2=zeros(1,length)

mu=zeros(1,length)

for i=1:length%i从1开始

mu(i)=(pi0*(p0.^x(i)))*((1-p0).^(1-x(i)))/((pi0*(p0.^x(i)))*((1-p0).^(1-x(i)))+((1-pi0)*(q0.^x(i)))*((1-q0).^(1-x(i))))

a(i)=mu(i)*x(i)

b(i)=(1-mu(i))*x(i)

end

pi=sum(mu)/length

p1=sum(a)/sum(mu)

p2=sum(b)/sum(1-mu)

这个是我刚开始学习EM算法时候写的，希望对你有帮助。

%%clear the way

% author : liuweimin

% ictcas

close all

clear

clc

%% settings

M=3 % number of Gaussian

N=1000% total number of data samples

th=1e-3 % convergent threshold

Nit=2000 % maximal iteration

Nrep=10% number of repetation to find global maximal

K=2 % demention of output signal

pi=3.141592653589793% in case it is overwriten by smae name variable

cond_num =1000 % prevent the singular covariance matrix in simulation data

plot_flag=1

print_flag=1

%% paramethers for random signal genrator

% random parameters for M Gaussian signals

mu_real = randn(K,M)% mean

cov_real =zeros(K,K,M) % covariance matrix

covd_real=zeros(K,K,M) % covariance matrix decomposition

for cm=1:M

while 1

covd_real(:,:,cm)=randn(K,K)

cov_real(:,:,cm)=covd_real(:,:,cm)*covd_real(:,:,cm)'

if cond(cov_real(:,:,cm))>cond_num

continue

else

break

end

end

end

% probablilty of a channel being selected

a_real = abs(randn(M,1))

a_real = a_real/sum(a_real) % normlize

if print_flag==1

a_real%类别概率的真值

mu_real%均值的真值

cov_real%协方差的真值

end

%% generate random sample of Gaussian vectors

%m=randdist(1,N,[1:M],a_real) % selector

rand_num_a=rand()%产生随机数

rand_num_b=rand()%产生随机数

while rand_num_a>=rand_num_b

rand_num_a=rand()%产生随机数

rand_num_b=rand()%产生随机数

end

x=randn(K,N)

for c=1:round(rand_num_a*N)

%sel=m(c)

x(:,c)=covd_real(:,:,1)*x(:,c)+mu_real(1)

end

for c=round(rand_num_a*N)+1:round(rand_num_b*N)

%sel=m(c)

x(:,c)=covd_real(:,:,2)*x(:,c)+mu_real(2)

end

for c=round(rand_num_b*N)+1:(N)

%sel=m(c)

x(:,c)=covd_real(:,:,3)*x(:,c)+mu_real(3)

end

%% EM Algorothm

% loop

f_best=-inf

for crep=1:Nrep

c=1

% initial values of parameters for EM

a=abs(randn(M,1)) % randomly generated

a=a/sum(a)% normlize, such that sum(a_EM)=1

mu=randn(K,M)

cov =zeros(K,K,M) % covariance matrix

covd=zeros(K,K,M) % covariance matrix decomposition

for cm=1:M

while 1

covd(:,:,cm)=randn(K,K)

cov(:,:,cm)=covd(:,:,cm)*covd(:,:,cm)'

if cond(cov(:,:,cm))>cond_num

continue

else

break

end

end

end

% iteration to find local maxima

break_flag=0

while 1

a_old= a

mu_old= mu

cov_old=cov

fprintf(1,'calculating probability pmx...\n')

pause(0)

% pmx(m,x|param)

pmx=zeros(M,N)

for cm=1:M

cov_cm=cov(:,:,cm)

if cond(cov_cm) >cond_num

break_flag=1

end

inv_cov_cm=inv(cov_cm)

mu_cm=mu(:,cm)

for cn=1:N

%p_cm=exp(-0.5*(x(:,cn)-mu_cm)'*inv_cov_cm*(x(:,cn)-mu_cm))

p_cm=a(cm,:)*exp(-0.5*(x(:,cn)-mu_cm)'*inv_cov_cm*(x(:,cn)-mu_cm))%这里加上a

pmx(cm,cn)=p_cm

end

pmx(cm,:)=pmx(cm,:)/sqrt(det(cov_cm))

end

pmx=pmx*(2*pi)^(-K/2)

fprintf(1,'calculating conditional probability, p...\n')

pause(0)

% conditional probability p(m|x,param) for estimated parameters

p=pmx./kron(ones(M,1),sum(pmx))

fprintf(1,'updating parametres\n')

pause(0)

a = 1/N*sum(p')'

mu = 1/N*x*p'*diag(1./a)

for cm=1:M

a_cm=a(cm)

mu_cm=mu(:,cm)

tmp=x-kron(ones(1,N),mu_cm)

cov(:,:,cm)=1/N*(kron(ones(K,1),p(cm,:)).*tmp)*tmp'*diag(1./a_cm)

end

t=max([norm(a_old(:)-a(:))/norm(a_old(:))

norm(mu_old(:)-mu(:))/norm(mu_old(:))

norm(cov_old(:)-cov(:))/norm(cov_old(:))])

if print_flag==1

fprintf('c=%04d: t=%f\n',c,t)

c=c+1

end

if t<th

break

end

if c>Nit

disp('reach maximal iteration')

break

end

if break_flag==1

disp('***** break on singular covariance matrix *****')

break

end

end

f=sum(log(sum(pmx.*kron(ones(1,N),a))))

if f>f_best

a_best=a

mu_best=mu

cov_best=cov

f_best=f

end

end

!echo 真实的rand_num_a

rand_num_a

!echo 真实的rand_num_b

rand_num_b-rand_num_a

!echo 真实的rand_num_c

1-rand_num_b

!echo 迭代出来的a_best

a_best

!echo 真实的均值矩阵

mu_real

!echo 迭代出来的均值矩阵

mu_best

for cs=1:M

!echo 真实的协方差矩阵

cov_real(:,:,cs)

!echo 迭代出来的协方差矩阵

cov_best(:,:,cs)

end

%% plot all

% for 2D (K=2) only

x1_vect=-1:0.02:1

x2_vect=-1:0.02:1

px=zeros(length(x1_vect), length(x2_vect))

for c1=1:length(x1_vect)

for c2=1:length(x2_vect)

for cm=1:3

cov_real_cm=cov_real(:,:,cm)

mu_real_cm=mu_real(:,cm)

a_real_cm=a_real(cm)

x_cm=[x1_vect(c1)

x2_vect(c2)]

pm=a_real_cm*(2*pi)^(-0.5*K)*det(cov_real_cm)^(-0.5)*exp(-0.5*x_cm'*inv(cov_real_cm)*x_cm)

px(c1,c2)=px(c1,c2)+pm

end

end

end

px_hat=zeros(length(x1_vect), length(x2_vect))

for c1=1:length(x1_vect)

for c2=1:length(x2_vect)

for cm=1:3

cov_cm=cov(:,:,cm)

mu_cm=mu(:,cm)

a_cm=a(cm)

x_cm=[x1_vect(c1)

x2_vect(c2)]

pm=a_cm*(2*pi)^(-0.5*K)*det(cov_cm)^(-0.5)*exp(-0.5*x_cm'*inv(cov_cm)*x_cm)

px_hat(c1,c2)=px_hat(c1,c2)+pm

end

end

end

figure(1)clfhold on

hold on

title('理论图')

xlabel('x')

ylabel('y')

mesh(x1_vect, x2_vect, px)

figure(2)clfhold on

mesh(x1_vect, x2_vect, px_hat)

title('统计图')

xlabel('x')

ylabel('y')

---