EM演算法做系統辨識

用EM演算法做系統辨識，問題描述：

采集了一批輸入輸出資料 ，但不確定各個樣本資料分別來自于兩個子模型中的哪一個：
模型1： y=k1x+b1+v，
模型2： y=k2x+b2+w，
其中v和w分別為服從均值為0的正態分布的白噪聲干擾項，試利用樣本資料，基于EM演算法對模型1和模型2的引數進行辨識，

關于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的資料，并生成白噪聲，

% 生成程序

% 白噪聲
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，

% 初始化引數
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，
繼續迭代，直至結束，

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

完整代碼如下：

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.m

function E=fun(a,x,y)
x=x(:); 
y=y(:); 
Y=a(1)*x+a(2);
E=y-Y; %M檔案結束