construct_bnet_hmm.m

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

function [bnet,evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,gausskwadnodes,...
    names,onames,order,inv_o]=...
   construct_bnet_hmm(D,S,M,S_item,data)
%bnet: bnet constructed from bnet software of murphy (struct array)
%evidence_nodes: struct(4): NOT YET IMPLEMENTED: ARE TREATED AS PARTIALLY
    %OBSERVED EVIDENCE NODES
    %evidence_nodes.nodenrs vector of (ordered) nodenrs of nodes that are observed for all cases
    %evidence_nodes.respcat nr of resp categories (not equal to nodesizes,
    %%evidence_nodes.nodesizes 
    %evidence_nodes.data cell array of corresponding columns in dataset  
%partial_evidence_nodes: struct(3): 
    %partial_evidence_nodes.nodenrs vector of (ordered) nodenrs of nodes that are partially observed 
    %partial_evidence_nodes.nodesizes number of resp categories (equals
    %nodesizes in the graph
    %partial_evidence_nodes.data cell array of of corresponding rows in dataset  
%terminal_merged_nodes: struct(5): 
    %terminal_merged_nodes.nodenrs vector of (ordered) nodenrs of terminal (partially) observed nodes 
    %terminal_merged_nodes.nodesizes:vector containing the state space of
    %the nodes, can always be set to vector of ones
    %terminal_merged_nodes.data cell array of corresponding sets of columns in dataset 
    %that are merged into single node 
    %terminal_merged_nodes.nrvars vector of numbers of nodes merged
    %into each (ordered) terminal node
    %terminal_merged_nodes.respcat cell array of nr of respcat of merged
    %nodes
%hid_nodes struct(2)
    %hid_nodes.nodenrs vector of (ordered) nodenrs of nodes that are
    %hidden
    %hid_nodes.nodesizes
%gausskwadnodes: struct(2)
    %gausskwadnodes.nodenrs (ordered) nrs continuous normal hidden nodes that are discretised
                 
    %gausskwadnodes.nodesizes

%the following output is not needed but can be useful during debugging
%names: names of the original order of the nodes
%onames: names for the reordered nodes (reordering is done during the
    %construction of the graphical model)
%order: vector: gives for the ordered node 1...end the corresponding nr of
    %the unordered nodes: onames=names(order)
%inv_o: vector: gives for the unordered node 1...end the corresponding nr of
    %the ordered nodes names=onames(inv_o)

%

%default values
evidence_nodes.nodenrs =[];
evidence_nodes.nodesizes =[];
evidence_nodes.data ={};  
evidence_nodes.respcat =[];
partial_evidence_nodes.nodenrs =[];
partial_evidence_nodes.nodesizes =[];
partial_evidence_nodes.data ={}; 
terminal_merged_nodes.nodenrs =[];
terminal_merged_nodes.nodesizes =[];
terminal_merged_nodes.data ={};
terminal_merged_nodes.nrvars =[];
terminal_merged_nodes.respcat ={};
hid_nodes.nodenrs=[];
hid_nodes.nodesizes=[];
gausskwadnodes.nodenrs=[];
gausskwadnodes.nodesizes=[];


disp('it is assumed that the data have a ');
disp(' separate row for each case, and');
disp('the columns are variables organised as follows:');
disp('first all items of first measurement occasion');
disp('then all items of second measurement occasion, etc');
disp('pres a key to continue');
pause;


%name the variables, easier to debug

names=cell(2*D,1);

for i=1:D
    names{i}=eval(['''timepoint_',int2str(i),'''']);
end

for i=1:D
   names{D+i}=eval(['''Y_',int2str(i),'''']);
end

nr_nodes=length(names);
%We can specify the topology using a cell array as follows, where each row specifies a connection between two named nodes: 

connections=cell(2*D-1,2);
for i=1:D-1
    connections{i,1}=eval(['''timepoint_',int2str(i),'''']);
    connections{i,2}=eval(['''timepoint_',int2str(i+1),'''']);
end
for i=1:D
       connections{D-1+i,1}=eval(['''timepoint_',int2str(i),'''']);
       connections{D-1+i,2}=eval(['''Y_',int2str(i),'''']);
end
    % convert this cell array to an adjacency matrix using the following function: 
[adj ,onames,order ] = franks_mk_adj_mat(connections, names,1);

inv_o=inv_order(order);
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%alternatively, specify the adjacency matrix directly 
%%%%%%%%%%%%%%%%%%%%%%%%%%%



hid_nodes.nodenrs=strmatch('timepoint_', onames) ;
hid_nodes.nodesizes=S*ones(1,length(hid_nodes.nodenrs));
obs_nodes=strmatch('Y', onames);

%define terminal_merged_nodes
    %specify which terminal observed nodes are merged 
    %for unordered nodes, each of the D Ys contains M observables
    unord_merged_nodes=cell(D,1);
    unord_merged_nodes_respcat=cell(D,1);
    unord_nr_merged_nodes=M*ones(D,1);
    for i=1:D
        unord_merged_nodes{i}=(i-1)*M+[1:M];
        unord_merged_nodes_respcat{i}=S_item*ones(M,1);
    end
terminal_merged_nodes.nodenrs=obs_nodes;

[y,ind]=sort(order(terminal_merged_nodes.nodenrs));
terminal_merged_nodes.respcat(ind)=unord_merged_nodes_respcat;
merged_nodes(ind)=unord_merged_nodes;
for i=1:D 
    terminal_merged_nodes.data{i}=(data(:,merged_nodes{i}))';
end

unord_nr_merged_nodes=M*ones(D,1);
nr_merged_nodes(ind)= unord_nr_merged_nodes; 
terminal_merged_nodes.nrvars=nr_merged_nodes;
terminal_merged_nodes.nodesizes=ones(D,1);



%define bnet
discrete_nodes = 1:nr_nodes;
node_sizes(terminal_merged_nodes.nodenrs) = 1;%for HMM:type models: either they are observed, and actual node size is then one
%or they are unobserved and thus can be left out because terminal hidden nodes have no influence
node_sizes(hid_nodes.nodenrs)=S;

%murphy's
bnet = mk_bnet(adj, node_sizes, 'discrete', discrete_nodes, 'names','onames');