EM算法作系统辨识

用EM算法作系统辨识，问题描述：javascript

采集了一批输入输出数据，但不肯定各个样本数据分别来自于两个子模型中的哪个：
模型1： y=k1x+b1+v，
模型2： y=k2x+b2+w，
其中v和w分别为服从均值为0的正态分布的白噪声干扰项。试利用样本数据，基于EM算法对模型1和模型2的参数进行辨识。

java

关于EM算法的理解能够看这篇文章硬币的例子https://blog.csdn.net/v_JULY_v/article/details/81708386算法

matlab源码见个人另外一篇，也可之间在下方代码复制。
https://download.csdn.net/download/weixin_42496224/13077074
1.数据生成
生成40%模型1和60%模型2的数据，并生成白噪声。

spa

% 生成过程

% 白噪声
x1 = randn(400,1);
x2 = randn(600,1);

% 数据生成

N = 1000;
x = zeros(N,1);

num_x1=1;
num_x2=1;

for i = 1 : N*0.4
    x(i) = i;
    y(i) = x(i)+1+x1(i);
end
for i = 1:0.6*N
    x(i+400) = i+400;
    y(i+400) = 2*x(i+400)+3+x2(i);
end

2.EM算法初始化
初始化中随意选取k1,b1,k2,b2
x1_para表示k1,b1;
x2_para表示k2,b2。

.net

% 初始化参数
x1_para = [1 2]';
x2_para = [3 3]';
x1_M_calulate = [];
x2_M_calulate = [];
y1_M_calulate = [];
y2_M_calulate = [];
M1_num = 1;
M2_num = 1;
% z表示x(i)的类别
z=[];

3.循环
循环分为E-step和M-step。
在E-step中，根据先前估计出的k1,b1,k2,b2分别计算出每一个点的y1和y2，比较|y-y1|和|y-y2|哪一个更小，小就表明当前点属于该模型的几率更大。
在M-step中，因为在前一步E-step中已经获得了每一个点更有可能属于的模型，将两个模型的全部点做非线性最小二乘拟合，拟合出新的k1,b1,k2,b2。
继续迭代，直至结束。

code

for o=1:100
% E-step
    M1_num=1;
    M2_num=1;
    clear x1_M_calulate;
    clear x2_M_calulate;
    clear y1_M_calulate;
    clear y2_M_calulate;
    x1_M_calulate = [];
    x2_M_calulate = [];
    y1_M_calulate = [];
    y2_M_calulate = [];
    for t=1:1000
        compare1 = abs(y(t)-x1_para(1)*x(t)-x1_para(2));
        compare2 = abs(y(t)-x2_para(1)*x(t)-x2_para(2));

        if compare1<compare2
            z(t)=1;
            x1_M_calulate(M1_num) = x(t);
            y1_M_calulate(M1_num) = y(t);
            M1_num = M1_num+1;
        else
            z(t)=2;
            x2_M_calulate(M2_num) = x(t);
            y2_M_calulate(M2_num) = y(t);
            M2_num = M2_num+1;
        end
    end

% M-step
    a0=[1 1]; 
    options=optimset('lsqnonlin'); 
    p1=lsqnonlin(@fun,a0,[],[],options,x1_M_calulate',y1_M_calulate');
    p2=lsqnonlin(@fun,a0,[],[],options,x2_M_calulate',y2_M_calulate');
    
    x1_para(1) = p1(1);
    x1_para(2) = p1(2);
    x2_para(1) = p2(1);
    x2_para(2) = p2(2);

end

完整代码以下：blog

clear;

clc;

% 设40%为y=x+1
% 60%为y=2x+3;
% 取1000个点;

%% 
% 生成过程

% 白噪声
x1 = randn(400,1);
x2 = randn(600,1);

% 数据生成

N = 1000;
x = zeros(N,1);

num_x1=1;
num_x2=1;

for i = 1 : N*0.4
    x(i) = i;
    y(i) = x(i)+1+x1(i);
end
for i = 1:0.6*N
    x(i+400) = i+400;
    y(i+400) = 2*x(i+400)+3+x2(i);
end

%%
% EM算法流程

% 初始化参数
x1_para = [1 2]';
x2_para = [3 3]';
x1_M_calulate = [];
x2_M_calulate = [];
y1_M_calulate = [];
y2_M_calulate = [];
M1_num = 1;
M2_num = 1;
% z表示x(i)的类别
z=[];

for o=1:100
% E-step
    M1_num=1;
    M2_num=1;
    clear x1_M_calulate;
    clear x2_M_calulate;
    clear y1_M_calulate;
    clear y2_M_calulate;
    x1_M_calulate = [];
    x2_M_calulate = [];
    y1_M_calulate = [];
    y2_M_calulate = [];
    for t=1:1000
        compare1 = abs(y(t)-x1_para(1)*x(t)-x1_para(2));
        compare2 = abs(y(t)-x2_para(1)*x(t)-x2_para(2));

        if compare1<compare2
            z(t)=1;
            x1_M_calulate(M1_num) = x(t);
            y1_M_calulate(M1_num) = y(t);
            M1_num = M1_num+1;
        else
            z(t)=2;
            x2_M_calulate(M2_num) = x(t);
            y2_M_calulate(M2_num) = y(t);
            M2_num = M2_num+1;
        end
    end

% M-step
    a0=[1 1]; 
    options=optimset('lsqnonlin'); 
    p1=lsqnonlin(@fun,a0,[],[],options,x1_M_calulate',y1_M_calulate');
    p2=lsqnonlin(@fun,a0,[],[],options,x2_M_calulate',y2_M_calulate');
    
    x1_para(1) = p1(1);
    x1_para(2) = p1(2);
    x2_para(1) = p2(1);
    x2_para(2) = p2(2);

end
x1_para
x2_para

另外建立一个文件命名fun.mip

function E=fun(a,x,y)
x=x(:); 
y=y(:); 
Y=a(1)*x+a(2);
E=y-Y; %M文件结束