V_n=MGraph_GaussEstimate_variance(sample_var,set_a,set_b)

function V_n=MGraph_GaussEstimate_variance(sample_var,set_a,set_b)
%Input: sample_var is the original sample variance matrix
%       set_a is the subsets {a1,a2,a3, ...am} of K={1 ,2,3 ,4...,k} , k is the number of vertix.
%       set_b is the complement of a in K, b=K\a;
%From book P/184, Graphical models in applied multiariate statistics
%This version is correct!!! MArch 26. 2002
%
%Example: sample_var=[10 1 5 4; 1 10 2 6; 5 2 10 3; 4 6 3 10]
%   set_a=[ 1 2; 2 3; 3 4 ; 1 4], set_b=[3 4; 1 4; 1 2; 2 3]
%K={1,2,3,4}
%Author: Junbai Wang

lnofset=size(set_a,1);
ii=1;
V_n=eye(size(sample_var)); %inital;
%added 09.2003
iscell_set=iscell(set_b);
%end added

while ii<5000
    old_Vn=V_n;
    for i=1:lnofset
        %added jbw 09.2003
        if ~iscell_set
            xa=unique(set_a(i,:));
            xb=unique(set_b(i,:));
        else
            xa=unique(set_a(i,:));;
            xb=unique(set_b{i});
        end
        %end added 2003
        %pause
   %partition of a and b sets
        var_xa=V_n(xa,xa); %sample_var(xa,xa);
        var_xb=V_n(xb,xb); %sample_var(xb,xb);
        cov_xbxa=V_n(xb,xa); %sample_var(xb,xa);
        cov_xaxb=V_n(xa,xb); %sample_var(xa,xb);
        B_ba=cov_xbxa/var_xa;
        V_bba=var_xb-cov_xbxa/var_xa*cov_xaxb;
        V_aa=sample_var(xa,xa); 
        V_aaUp=V_aa*B_ba';
        V_aaLow=B_ba*V_aa;
        V_aaOther=V_bba+B_ba*V_aa*B_ba';

        V_n(xa,xa)=V_aa;
        V_n(xa,xb)=V_aaUp;
        V_n(xb,xa)=V_aaLow;
        V_n(xb,xb)=V_aaOther   ;
   end %i
   max_device=max(max(abs(old_Vn-V_n)));
   if max_device<0.0001
       break;
   end
   ii=ii+1;
end