www.pudn.com > RPEM_Source_Code.zip > rpem_adaptive.m

% Adaptive RPEM algorithm 
 
function [new_betta, new_m, new_S, betta_trace, alpha_trace, m_trace, S_trace] = rpem_adaptive(x, k, lrate, betta, m, S) 
total_epoch_no = input('Please input the number of epoches [default: 10]'); 
if isempty(total_epoch_no) 
        total_epoch_no = 10; 
end 
epoch_no = total_epoch_no; 
 
[d, N] = size(x); 
 
% Parameter  
h = zeros(1, k); 
alpha = zeros(1, k); 
px = zeros(1, k); 
new_betta = zeros(1, k); 
new_m = zeros(d, k); 
new_S = zeros(d,d,k); 
z = zeros(d, k); 
 
flag = 1; 
epoch_count = 0; 
colorset_m = zeros(d*k, 3); 
 
mf = figure; 
if d == 2 
    xf = figure; 
end 
 
 % Record the parameters 
 alpha_trace(:,epoch_count+1) = vec(exp(betta) / sum(exp(betta))); 
 betta_trace(:, epoch_count+1) = vec(betta); 
 m_trace(:,epoch_count+1) = vec(m); 
 S_trace(:,epoch_count+1) = vec(S); 
  
% Set the learning rate of betta 
lrate_betta = 0.1 * lrate; 
 
while (flag == 1) 
    for t = 1:N 
        % Perform Step 1 
        max_betta = max(betta); 
        if max_betta > 100 
           betta = betta - 100; 
       end 
        alpha = exp(betta)/sum(exp(betta)); 
 
 
        for j = 1:k 
            z(:,j) = x(:,t)-m(:,j); 
            px(j) = exp(-0.5*z(:,j)'*S(:,:,j)*z(:,j))*sqrt(det(S(:,:,j))); 
        end 
        h = (alpha.*px)/sum(alpha.*px); 
        [max_h, c] = max(h); 
         
        % Perform Step 2 
        for j = 1:k 
                if j == c 
                    temp_g = 2 - h(c); 
                else 
                    temp_g = -h(j); 
                end 
 
                if j == c 
                          new_betta(j) = betta(j) + lrate_betta*(2-h(c)-alpha(c));                             
                else 
                          new_betta(j) = betta(j) - lrate_betta*(h(j)+alpha(j));                             
                end                 
                new_m(:,j) = m(:,j) +  lrate*temp_g*S(:,:,j)*z(:,j); 
			    new_S(:,:,j) = (1+0.5*lrate*temp_g)*S(:,:,j) - 0.5*lrate*temp_g*S(:,:,j)*z(:,j)*z(:,j)'*S(:,:,j); 
       end   
       betta = new_betta; 
       m = new_m; 
       S = new_S; 
  end 
 epoch_no = epoch_no - 1; 
 epoch_count = epoch_count + 1; 
 fprintf('Epoch = %d \n', epoch_count); 
 
 % Record the parameters 
 alpha_trace(:,epoch_count+1) = vec(exp(betta) / sum(exp(betta))); 
 betta_trace(:, epoch_count+1) = vec(betta); 
 m_trace(:,epoch_count+1) = vec(m); 
 S_trace(:,epoch_count+1) = vec(S); 
  
  if epoch_no == 0 
            disp('The learned alpha are:'); 
            alpha = exp(new_betta)/sum(exp(new_betta)) 
            disp('The learned positions of m_js are:'); 
             new_m 
            disp('The learned sigma are:'); 
            for j = 1:k 
               inv(new_S(:,:,j)) 
            end 
                  
        	figure(mf); 
        	hold off; 
        	for t=1:d*k 
        		ploth = plot(0:epoch_count, m_trace(t,:)); 
        		set(ploth, 'Color', colorset_m(t,:)); 
        		hold on; 
        	end 
        	title('Learning Curve of Parameter m_js'); 
        	xlabel('No. of data points scanned (x 10^3)'); 
        	drawnow; 
                   
            if d == 2 
        	    figure(xf); 
        	    hold off; 
        	    plot(x(1,:), x(2,:), '.'); 
                hold on; 
                plot(m(1,:), m(2,:), 'black*'); 
        	    title('Positions of Parameter m_js in Input Space'); 
        	    drawnow; 
            end 
            flag_c = input('Do you want to continue to learning the parameters (yes:1, no:0) [default: 1]?'); 
            if isempty(flag_c) 
                    flag_c = 1; 
           end 
           if flag_c ~= 1 
               flag = 0; 
           else 
               epoch_no = input('Please input the number of epoches [default 10]'); 
               if isempty(epoch_no) 
                    epoch_no = 10; 
               end 
               total_epoch_no = total_epoch_no + epoch_no; 
           end 
       end 
   end