Daniel@0: function seq = sample_dbn(bnet, varargin)
Daniel@0: % SAMPLE_DBN Generate a random sequence from a DBN.
Daniel@0: % seq = sample_dbn(bnet, ...)
Daniel@0: %
Daniel@0: % seq{i,t} contains the values of the i'th node in the t'th slice.
Daniel@0: %
Daniel@0: % Optional arguments:
Daniel@0: %
Daniel@0: % length - length of sequence  to be generated (can also just use sample_dbn(bnet,T))
Daniel@0: % stop_test - name of a function which is used to decide when to stop;
Daniel@0: %   This will be called as feval(stop_test, seq(:,t))
Daniel@0: %   i.e., stop_test is passed a cell array containing all the nodes in the current slice.   
Daniel@0: % evidence - initial evidence; if evidence{i,t} is non-empty, this node won't be sampled.
Daniel@0: 
Daniel@0: args = varargin;
Daniel@0: nargs = length(args);
Daniel@0: 
Daniel@0: if (nargs == 1) & ~isstr(args{1})
Daniel@0:   % Old syntax: sample_dbn(bnet, T)
Daniel@0:   T = args{1};
Daniel@0: else
Daniel@0:   % get length
Daniel@0:   T = 1;
Daniel@0:   for i=1:2:nargs
Daniel@0:     switch args{i},
Daniel@0:      case 'length',      T = args{i+1}; 
Daniel@0:      case 'evidence',    T = size(args{i+1}, 2);
Daniel@0:     end
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: ss = length(bnet.intra);
Daniel@0: % set default arguments
Daniel@0: seq = cell(ss, T);
Daniel@0: stop_test = [];
Daniel@0: for i=1:2:nargs
Daniel@0:   switch args{i},
Daniel@0:    case 'evidence',    seq = args{i+1}; % initialise observed nodes
Daniel@0:    case 'stop_test',  stop_test = args{i+1};
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: t = 1;
Daniel@0: for i=1:ss
Daniel@0:   if ~isempty(stop_test) | isempty(seq{i,t})
Daniel@0:     ps = parents(bnet.dag, i);
Daniel@0:     e = bnet.equiv_class(i,1);
Daniel@0:     pvals = seq(ps);
Daniel@0:     seq{i,t} = sample_node(bnet.CPD{e}, pvals);
Daniel@0:     %fprintf('sample i=%d,t=%d,val=%d,ps\n', i, t, seq(i,t)); pvals(:)'
Daniel@0:   end
Daniel@0: end
Daniel@0: t = 2;
Daniel@0: done = 0;
Daniel@0: while ~done
Daniel@0:   for i=1:ss
Daniel@0:     if ~isempty(stop_test) | isempty(seq{i,t})
Daniel@0:       ps = parents(bnet.dag, i+ss) + (t-2)*ss;
Daniel@0:       e = bnet.equiv_class(i,2);
Daniel@0:       pvals = seq(ps);
Daniel@0:       seq{i,t} = sample_node(bnet.CPD{e}, pvals);
Daniel@0:       %fprintf('sample i=%d,t=%d,val=%d,ps\n', i, t, seq(i,t)); pvals(:)'
Daniel@0:     end
Daniel@0:   end
Daniel@0:   if ~isempty(stop_test)
Daniel@0:     done = feval(stop_test, seq(:,t));
Daniel@0:   else
Daniel@0:     if t==T
Daniel@0:       done = 1;
Daniel@0:     end
Daniel@0:   end
Daniel@0:   t = t + 1;
Daniel@0: end