num_infomatrix_anal_score.m

The BNL toolbox is a set of Matlab functions for defining and estimating the parameters of a Bayesi

function info=num_infomatrix_anal_score(link,parms,restparms,design,parents,node_sizes,equiv_class,...
    postorder,septable, cliquetable,preorder,equiv_class_time,clq_ass_to_node,evidence_nodes,...
    partial_evidence_nodes,terminal_merged_nodes,hid_nodes,gausskwadnodes,N,e)
%computes infomatrix using finite differences
%score function is derived analytically (score for complete data = score
%for incomplete data)
if nargin==19 e=.0001; end
info=[];
f_d2=[];

pot_to_CPT=link_pot_to_CPT(cliquetable,clq_ass_to_node,parents);
disp('computing info matrix by numerical differentation of score function');

for i=1:length(parms)
   i
    p=parms{i};
    f=find(restparms{i}~=1);
    if ~isempty(f);
    for j=1:length(f)
        parms2=parms;
        parms3=parms;
        parms2{i}(f(j))=parms{i}(f(j))+e;
        parms3{i}(f(j))=parms{i}(f(j))-e;
        i;
        j;
        lin_pred=construct_lin_pred(parms2,design,parents,node_sizes,equiv_class,equiv_class_time,terminal_merged_nodes,N);
        equiv_class_CPTs=construct_equiv_class_CPT(link,lin_pred,N);
        bigCPTs=construct_bigCPTs(equiv_class_CPTs,equiv_class,...
        evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,...
        gausskwadnodes,N);
        [cliquetable, septable,schaal_init]=franks_init_pot(cliquetable,septable,bigCPTs,parents);
        %collect
        [septable,cliquetable,schaal_collect]=propagate_messages(postorder,septable, cliquetable);%kan sneller!
        %distribute
        [septable,cliquetable,schaal_distribute]=propagate_messages(preorder,septable, cliquetable);
        joint_prob_tabs=compute_JPTs(cliquetable,pot_to_CPT,parents,node_sizes,N);
        %compute score
        scor2=score(link,parms2,restparms,joint_prob_tabs,design,equiv_class,equiv_class_time,terminal_merged_nodes,N);
        
        
        lin_pred=construct_lin_pred(parms3,design,parents,node_sizes,equiv_class,equiv_class_time,terminal_merged_nodes,N);
        equiv_class_CPTs=construct_equiv_class_CPT(link,lin_pred,N);
        bigCPTs=construct_bigCPTs(equiv_class_CPTs,equiv_class,...
        evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,...
        gausskwadnodes,N);
        [cliquetable, septable,schaal_init]=franks_init_pot(cliquetable,septable,bigCPTs,parents);
        %collect
        [septable,cliquetable,schaal_collect]=propagate_messages(postorder,septable, cliquetable);%kan sneller!
        %distribute
        [septable,cliquetable,schaal_distribute]=propagate_messages(preorder,septable, cliquetable);
        joint_prob_tabs=compute_JPTs(cliquetable,pot_to_CPT,parents,node_sizes,N);
        %compute score
        scor3=score(link,parms3,restparms,joint_prob_tabs,design,equiv_class,equiv_class_time,terminal_merged_nodes,N);
        
        f_d2=[f_d2 ((scor2-scor3)/(2*e))];
    end
    end
end


info=-f_d2;