%% 程序说明
%% 该程序为混合正态分布参数估计算法
%% 微博：iyx_yao
%% 欢迎交流

%% 产生拟合数据
k = 3; % 高斯分布的个数
mu = [0  50  100];
sigma = [10 10 10];
nk = [1/3 1/3 1/3];
x = [];
for i = 1:k
    x = [x normrnd(mu(i),sigma(i),1,300)];
end


%% 参数初始化：kmeans算法估计
idx = kmeans(x',k);
for i = 1:k
    N(i) = length(x(idx == i));
    a(i) = N(i)/ length(x);
    u(i) = mean(x(idx == i));
    o(i) = std(x(idx == i));
end

%% 迭代程序
t = 1;
while t < 100
    Es =0;
    for i =1:k
    Es = Es +a(i) * normpdf(x,u(i),o(i));
    end
    for i =1:k
    E = a(i) *normpdf(x,u(i),o(i));
    N(i) = sum(E./Es);
    a(i) = N(i)/length(x);
    u(i) =  sum(E./Es.*x)/N(i);
    o(i) = sqrt(sum(E./Es.*(x-u(i)).^2)/N(i));
    end
    t = t +1;
end
disp('理论值')
[nk,mu,sigma]
disp('计算值')
[a,u,o]