| [old_parcorrf,delet_edge,record_edge]=MGraph_GaussModel(old_var,N,q,cutoff,deviceChoice)
|
function [old_parcorrf,delet_edge,record_edge]=MGraph_GaussModel(old_var,N,q,cutoff,deviceChoice)
global graph_form
%%Graphical GAussian model for continues data
%Input:
%old_var, covariance matrix of inital data;
%N, number of observations
%q, number of variables = size of var
%cutoff, significant level 0.05 defaule
%
%this program is seems ok !! MArch 26.2002
%but more test run need do
%Author Junbai Wang
t0=clock;
%%test march 20 2002
%P/181, Graphical models in applied multivariate statistics
%old_var=[ 302.29 125.78 100.43 105.07 116.07; ...
% 125.78 170.88 84.19 93.6 97.89;...
% 100.43 84.19 111.6 110.84 120.49; ...
% 105.07 93.6 110.84 217.88 153.77; ...
% 116.07 97.89 120.49 153.77 294.37]
%N=88 % number of observations
%q=5 % number of variables
%cutoff=0.05;
sigma_X=[];
delet_edge=1;
record_edge=zeros(1,2);
oldest_var=old_var;
[nrow ncol]=size(old_var);
totalEdges=sum(sum(triu(ones(nrow,nrow))))-nrow;
%Make a completed set of a initialize of set_a only up part of the matris
ii=1;
for i=1:nrow
for j=i+1:ncol
set_a(ii,:)=[i,j];
ii=ii+1;
end
end
set_K=1:nrow; %number of vertix
temp_parcorrf=ones([nrow ncol]);
%make the full model
oldModel=MGraph_num2string(1:q);
temp_newModel=oldModel;
while 1
%start to fit the models and find the minume partial corrcoef in initial modle
%until the device detection is full then stop
old_corrf=old_var./sqrt(diag(old_var)*diag(old_var)');
inv_old_var=pinv(old_var); %bug at here sometime the inverse of variance equal inf Dec 22. 2002
old_parcorrf=-inv_old_var./sqrt(abs(diag(inv_old_var)*diag(inv_old_var)'));
if delet_edge==1
oldest_var=old_var;
end
%find min parcorrf
min_par=sort(abs(old_parcorrf(find(abs(old_parcorrf)>0.000001))));
for ki=1:length(min_par)
[min_pari, min_parj]=find(abs(old_parcorrf)==min_par(ki));
%add 12.03 for check decomposable modle
%make new subgroup
if graph_form==1
delet_edge_location=[min_pari,min_parj];
delet_edge_string=MGraph_num2string(delet_edge_location);
new_model=MGraph_makeSubGraph(temp_newModel,delet_edge_string);
clear numberVector
for new_idx=1:length(new_model)
numberVector{new_idx}=num2str(MGraph_string2num(new_model(new_idx)));
end
input_k_cliques=(numberVector)';
end
if graph_form==1 %decomposable modle
is_SDordering=MGraph_isSDordering(input_k_cliques);
elseif graph_form==2 %unrestricted model
is_SDordering=1;
else
error('Please chose type of graphical model in ''MGraph properties''! Program Stop!')
break;
end
%end added
if is_SDordering==1 %if it is a decomposiable model
if (~ismember([ min_pari min_parj ],record_edge,'rows') & ~ismember([ min_parj min_pari ],record_edge,'rows'))
min_pari=min_pari(1);
min_parj=min_parj(1);
break;
end
end
end
%replace min parcorrf as zero, at here , we estimate the fitted covariance in this point
%But I have to check how to transfer updated parcorrcoef to corrcoef or covariance matrix!!
set_b=[min_pari,min_parj] ;
size(set_a), size(set_b);
set_a=setdiff(set_a,set_b,'rows'); % new set of a excluded the min pair of parcorrelations
set_b=[min_parj,min_pari] ;
set_a=setdiff(set_a,set_b,'rows'); % new
%make complement set of set_bb to the new set_a
[set_row set_col]=size(set_a) ; % set_Col should always be 2
for ki=1:set_row
set_bb(ki,:)=setdiff(set_K,set_a(ki,:));
end
fit_var=MGraph_GaussEstimate_variance(old_var,set_a,set_bb);
inv_fit_var=pinv(fit_var);
fit_parcorrf=-inv_fit_var./sqrt(diag(inv_fit_var)*diag(inv_fit_var)');
%there is a bug here nov 14.2002
%the choose of device change the final results, the old_var device is better
%d=abs(N*log(det(old_parcorrf)/det(fit_parcorrf)))
%N
if deviceChoice==1
det_fit=det(fit_var)
det_old=det(old_var)
%det_old=det(old_parcorrf)
%det_fit=det(fit_parcorrf)
%if (det_fit<0.0001 & det_old<0.0001 )
% d=1;
%else
log(det_old),log(det_fit)
d=abs(N*(log(det_old)-log(det_fit)));
%end
elseif deviceChoice==2
%%new device added at June 2003 %cutoff is 0.95 to 0.96
fit_corrf=fit_var./sqrt(diag(fit_var)*diag(fit_var)');
dX0=1-old_corrf;
delt=1-fit_corrf;
d=(sum(sum(triu(delt.^2))+sum(triu(dX0.^2))-2*sum(triu(delt).*triu(dX0)))./sum(sum(triu(delt.^2))));
elseif deviceChoice==3
%F-test
det_fit=det(fit_var)
det_old=det(oldest_var)
d=abs((det_fit/det_old));
elseif deviceChoice==4
%book Mathematical statistics with applicaitons, Mendenhall, Page 443
det_fit=det(fit_var)
det_old=det(oldest_var)
d=abs(N-delet_edge-1)*(det_fit/det_old);
elseif deviceChoice==5
det_fit=det(fit_var)
det_old=det(old_var)
d=abs(N*(log(det_old)-log(det_fit)));
k=totalEdges-delet_edge;
f=(N-k)*(exp(d/N)-1);
elseif deviceChoice==6
%Bartlett correction
det_fit=det(fit_var)
det_old=det(old_var)
d=abs(N*(log(det_old)-log(det_fit)));
d=d*(1-1/N*(N-totalEdges+1+1/2*(1+1+1)));
end
marray_debuge(['Device= ',num2str(d)]);
if deviceChoice==3
df1=N-delet_edge-1
df2=N-1
df=df1;
elseif deviceChoice==4
df=N-1-delet_edge;
elseif deviceChoice==5
df1=1;
df2=N-k;
df=df1;
else
df=1 ;%degree of freedom
end
marray_debuge(['Degree of freedom= ', num2str(df)]);
sigma_X(delet_edge)=d;
if deviceChoice==3
chtest(delet_edge)=marray_FishF(d,df1,df2);
elseif deviceChoice==4
chtest(delet_edge)=Chisq(d,df);
elseif deviceChoice==5
chtest(delet_edge)=marray_FishF(f,df1,df2);
else
chtest(delet_edge)=Chisq(d,df);
end
marray_debuge(['Chi square test P= ',num2str(chtest)]);
marray_debuge(['Device = ',num2str(sigma_X)]);
%pause
if d==0 & chtest(delet_edge)==1
old_parcorrf=old_parcorrf; %.*temp_parcorrf
delet_edge=delet_edge-1;
break;
end
if chtest(delet_edge) <cutoff
old_parcorrf=old_parcorrf; %.*temp_parcorrf
delet_edge=delet_edge-1;
break;
else
old_var=fit_var;
%old_parcorrf=fit_parcorrf; %updated covancice and delete the edge
record_edge(delet_edge,:)=[min_pari,min_parj];
if graph_form==1
temp_newModel=new_model ;
end
marray_debuge('Selected model:');
marray_debuge(strvcat(temp_newModel));
marray_debuge(['Deleted edge ']);
marray_debuge([num2str(record_edge(end,:))]);
marray_debuge(' ');
temp_parcorrf(min_pari,min_parj)=0;
temp_parcorrf(min_parj,min_pari)=0;
delet_edge=delet_edge+1;
end
end
%temp_newModel
etime(clock,t0)
|
|
