www.pudn.com > RVM1.zip > SeqSparBayeLear.m, change:2011-06-22,size:2815b


% ¦µ = [¦Õ1 . . . ¦ÕM ] is the N ¡Á M ¡®design¡¯ matrix whose 
% columns comprise the complete set of M ¡®basis vectors¡¯. 
function [mu model_idx] = SeqSparBayeLear(basis,target,maxIter) 
[N M] = size(basis); 
infi = 10^6; 
alpha = ones(M,1)*infi; 
model_idx = zeros(1,M); 
%% innitialise sigma_s to var[target]*0.1 
sigma_s = var(target)*0.1; 
%% initialise alpha with a single basis vector  
alpha(1) = norm(basis(:,1))^2/(norm((basis(:,1)'*target))^2/norm(basis(:,1))^2-sigma_s); 
model_idx(1) = 1; 
%% caculate ¦²  
idx = find(model_idx == 1);  
A = diag(alpha(idx)); 
su_sigma = (A+basis(:,idx)'*basis(:,idx)/sigma_s)^(-1); 
%% caculate  u 
mu = su_sigma*basis(:,idx)'*target/sigma_s; 
%% initialise s q 
s = zeros(M,1); 
q = zeros(M,1); 
B = 1/sigma_s*eye(N,N); 
MA = B*basis(:,idx)*su_sigma*basis(:,idx)'*B; 
for m = 1:M 
    Sm = basis(:,m)'*B*basis(:,m)-basis(:,m)'*MA*basis(:,m); 
    Qm = basis(:,m)'*B*target-basis(:,m)'*MA*target; 
    if alpha(m) < infi 
    s(m) = alpha(m)*Sm/(alpha(m)-Sm); 
    q(m) = alpha(m)*Qm/(alpha(m)-Sm); 
    else 
        s(m) = Sm; 
        q(m) = Qm; 
    end 
end 
iter = maxIter; 
while iter ~=0 
    iter = iter -1 
    %% compute  ¦Èi  
    i = mod(iter,M)+1;%randi(M); 
    sitai = q(i)^2-s(i); 
    %% model  
    old_idx = find(model_idx==1); 
    old_basis = basis(:,old_idx); 
    old_alpha = alpha(old_idx); 
    if sitai>0 && alpha(i)<infi 
        %re-estimate alpha(i) 
        alpha(i) = s(i)^2/(q(i)^2-s(i)); 
    end 
    if sitai>0 && alpha(i)==infi 
        model_idx(i) = 1; 
        % update alpha(i) 
        alpha(i) = s(i)^2/(q(i)^2-s(i)); 
    end 
    if sitai<=0 && alpha(i)<infi 
        if model_idx(i) ~=0 
        model_idx(i) = 0; 
        alpha(i) = infi; 
        end 
    end 
    %% updata ¦Ò2  
    sigma_s = norm(target-old_basis*mu)^2/(N-M+old_alpha'*diag(su_sigma)); 
    %% update ¦² ? 
    idx = find(model_idx==1); 
    A = diag(alpha(idx)); 
    su_sigma = (A+basis(:,idx)'*basis(:,idx)/sigma_s)^(-1); 
    mu = su_sigma*basis(:,idx)'*target/sigma_s; 
    %% updata sm and qm 
    B = 1/sigma_s*eye(N,N); 
    MA = B*basis(:,idx)*su_sigma*basis(:,idx)'*B; 
    for k = 1:M 
        Sm = basis(:,m)'*B*basis(:,m)-basis(:,m)'*MA*basis(:,m); 
        Qm = basis(:,m)'*B*target-basis(:,m)'*MA*target; 
        if alpha(m) < infi 
            s(m) = alpha(m)*Sm/(alpha(m)-Sm); 
            q(m) = alpha(m)*Qm/(alpha(m)-Sm); 
        else 
            s(m) = Sm; 
            q(m) = Qm; 
        end 
    end 
    %% conergence 
        sita = q.^2-s; 
 if max(sita(model_idx == 0))<=0 
    if size(old_idx) == size(idx) 
        if old_idx == idx 
            change_alpha = abs(log(old_alpha)-log(alpha(idx))); 
            if max(change_alpha) <= 10^-3; 
                break; 
            end 
        end 
    end 
 end 
end