rpem_main.m

EM算法

% ************************************************************************
% This algorithm is to perform Rival Penalized Expectation-Maximization
% (RPEM) Algorithm.
% Version 1 (November 12, 2002)
% ************************************************************************

clear all;
close all;

% Determine the number of observations
onum = input('Please input the number of observations [default: 1000]');
if isempty(onum)
	onum = 1000;
end

% Setting for generating 2-dimension data clusters
disp('1. Generate three 2-dimension ball-shaped clusters well separated');
disp('2. Generate three 2-dimension clusters well separated, but each cluster is in general shaped.');
disp('3. Generate three 2-dimension ball-shaped clusters that are overlapped moderately');
disp('4. Generate three 2-dimension clusters that are overlapped moderately, but each cluster is in general shaped.');
disp('5. Generate three 2-dimension clusters that are seriously overlapped, and each cluster is in general shaped.')
choice = input('Please input the choice (default: 1):');
if isempty(choice)
	choice = 1;
end

x_dim = 2;

cluster_no = 3;
range = [0.3, 0.7, 1];

mu = zeros(x_dim, cluster_no);
sigma = zeros(x_dim, x_dim, cluster_no);

if choice == 1
	mu = [1.0, 1.0, 5.0; 
      	  1.0, 5.0, 5.0];
	for j = 1:cluster_no
		sigma(:,:,j) = 0.1*eye(x_dim);
	end
elseif choice == 2
	mu = [1.0, 1.0, 5.0; 
           	1.0, 5.0, 5.0];
	sigma(:,:,1) = [0.1, 0.05;
                   0.05, 0.2];
    sigma(:,:,2) = 0.1*eye(x_dim);
    sigma(:,:,3) = [0.1, -0.05;
                  -0.05, 0.1];
elseif choice == 3
	mu = [1.0, 1.0, 2.5; 
      	1.0, 2.5, 2.5];
	for j = 1:cluster_no
		sigma(:,:,j) = 0.15*eye(x_dim);
	end
elseif choice == 4
	mu = [1.0, 1.0, 2.5; 
      	1.0, 2.5, 2.5];
	sigma(:,:,1) = [0.15, 0.05;
                   0.05, 0.2];
    sigma(:,:,2) = 0.15*eye(x_dim);
    sigma(:,:,3) = [0.15, -0.1;
                    -0.1, 0.15];
else
	mu = [1.0, 1.0, 2.5; 
      	1.0, 2.5, 2.5];
	sigma(:,:,1) = [0.2, 0.05;
                   0.05, 0.3];
    sigma(:,:,2) = 0.2*eye(x_dim);
    sigma(:,:,3) = [0.2, -0.1;
                   -0.1, 0.2];
end

% Generate data
x = gm_generation(x_dim, onum, range, mu, sigma);

% perform RPEM algorithm
k = input('Please input the number of seed points [default: k = 3]');
if isempty(k)
        k = 3;
end

% learning rate
lrate = 0.001;

% Initialize parameters
betta = zeros(1, k);
betta_init = betta;
m = zeros(x_dim, k);
S = zeros(x_dim, x_dim, k);
for j = 1:k
   	m(:,j) = x(:,j)+rand(x_dim,1);
    S(:,:,j) = eye(x_dim);
end
m_init = m;
S_init = S;
for j =1:k
        inv_S(:,:,j) = inv(S(:,:,j));
        inv_S_init = inv_S;
end


% Display the true parameters
disp('The true parameters are:');
disp('alpha = ');
[0.3, 0.4, 0.3]
disp('mu =')
mu
disp('sigma =')
sigma

if x_dim == 2
        xf = figure;
   	    figure(xf);
  	    hold off;
   	    plot(x(1,:), x(2,:), '.');
        hold on;
        plot(m_init(1,:), m_init(2,:), 'black*');
  	    title('Initial Positions of Parameter m_js in Input Space');
  	    drawnow;
end

% Learning algorithm 
[betta, m, inv_S, b_trace, alpha_trace, m_trace, S_trace] = rpem_adaptive(x, k, lrate, betta, m, inv_S);