wolffd@0: function bnet = mk_hhmm3(varargin)
wolffd@0: % MK_HHMM3 Make a 3 level Hierarchical HMM
wolffd@0: % bnet = mk_hhmm3(...)
wolffd@0: %
wolffd@0: % 3-layer hierarchical HMM where level 1 only connects to level 2, not 3 or obs.
wolffd@0: % This enforces sub-models (which differ only in their Q1 index) to be shared.
wolffd@0: % Also, we enforce the fact that each model always starts in its initial state
wolffd@0: % and only finishes in its final state. However, the prob. of finishing (as opposed to
wolffd@0: % self-transitioning to the final state) can be learned.
wolffd@0: % The fact that we always finish from the same state means we do not need to condition
wolffd@0: % F(i) on Q(i-1), since finishing prob is indep of calling context.
wolffd@0: %
wolffd@0: % The DBN is the same as Fig 10 in my tech report.
wolffd@0: %
wolffd@0: %   Q1 ---------->  Q1
wolffd@0: %   |              / |
wolffd@0: %   |             /  |
wolffd@0: %   |  F2 -------    |
wolffd@0: %   |  ^         \   |
wolffd@0: %   | /|          \  |
wolffd@0: %   v  |           v v
wolffd@0: %   Q2-| -------->   Q2
wolffd@0: %  /|  |             ^
wolffd@0: % / |  |            /|
wolffd@0: % | |  F3 ---------/ |
wolffd@0: % | |  ^           \ |
wolffd@0: % | v /              v
wolffd@0: % | Q3 ----------->  Q3
wolffd@0: % |  |    
wolffd@0: % \  | 
wolffd@0: %  v v    
wolffd@0: %   O
wolffd@0: %
wolffd@0: %
wolffd@0: % Optional arguments in name/value format [default]
wolffd@0: %
wolffd@0: % Qsizes      - sizes at each level [ none ]
wolffd@0: % Osize       - size of O node [ none ]
wolffd@0: % discrete_obs - 1 means O is tabular_CPD, 0 means O is gaussian_CPD [0]
wolffd@0: % Oargs       - cell array of args to pass to the O CPD  [ {} ]
wolffd@0: % transprob1  - transprob1(i,j) = P(Q1(t)=j|Q1(t-1)=i)  ['ergodic']
wolffd@0: % startprob1  - startprob1(j) = P(Q1(t)=j)  ['leftstart']
wolffd@0: % transprob2  - transprob2(i,k,j) = P(Q2(t)=j|Q2(t-1)=i,Q1(t)=k)  ['leftright']
wolffd@0: % startprob2  - startprob2(k,j) = P(Q2(t)=j|Q1(t)=k)  ['leftstart']
wolffd@0: % termprob2   - termprob2(j,f) = P(F2(t)=f|Q2(t)=j)  ['rightstop']
wolffd@0: % transprob3  - transprob3(i,k,j) = P(Q3(t)=j|Q3(t-1)=i,Q2(t)=k)  ['leftright']
wolffd@0: % startprob3  - startprob3(k,j) = P(Q3(t)=j|Q2(t)=k)  ['leftstart']
wolffd@0: % termprob3   - termprob3(j,f) = P(F3(t)=f|Q3(t)=j)  ['rightstop']
wolffd@0: %
wolffd@0: % leftstart means the model always starts in state 1.
wolffd@0: % rightstop means the model always finished in its last state (Qsize(d)).
wolffd@0: %
wolffd@0: % Q1:Q3 in slice 1 are of type tabular_CPD
wolffd@0: % Q1:Q3 in slice 2 are of type hhmmQ_CPD.
wolffd@0: % F2 is of type hhmmF_CPD, F3 is of type tabular_CPD.
wolffd@0: 
wolffd@0: ss = 6; D = 3;
wolffd@0: Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; obs = 6;
wolffd@0: Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
wolffd@0: names = {'Q1', 'Q2', 'Q3', 'F3', 'F2', 'obs'};
wolffd@0: 
wolffd@0: intra = zeros(ss);
wolffd@0: intra(Q1, Q2) = 1;
wolffd@0: intra(Q2, [F2 Q3 obs]) = 1;
wolffd@0: intra(Q3, [F3 obs]) = 1;
wolffd@0: intra(F3, F2) = 1;
wolffd@0: 
wolffd@0: inter = zeros(ss);
wolffd@0: inter(Q1,Q1) = 1;
wolffd@0: inter(Q2,Q2) = 1;
wolffd@0: inter(Q3,Q3) = 1;
wolffd@0: inter(F2,[Q1 Q2]) = 1;
wolffd@0: inter(F3,[Q2 Q3]) = 1;
wolffd@0: 
wolffd@0: 
wolffd@0: % get sizes of nodes
wolffd@0: args = varargin;
wolffd@0: nargs = length(args);
wolffd@0: Qsizes = [];
wolffd@0: Osize = 0;
wolffd@0: for i=1:2:nargs
wolffd@0:   switch args{i},
wolffd@0:    case 'Qsizes', Qsizes = args{i+1}; 
wolffd@0:    case 'Osize', Osize = args{i+1}; 
wolffd@0:   end
wolffd@0: end
wolffd@0: if isempty(Qsizes), error('must specify Qsizes'); end
wolffd@0: if Osize==0, error('must specify Osize'); end
wolffd@0:   
wolffd@0: % set default params
wolffd@0: discrete_obs = 0;
wolffd@0: Oargs = {};
wolffd@0: startprob1 = 'ergodic';
wolffd@0: startprob2 = 'leftstart';
wolffd@0: startprob3 = 'leftstart';
wolffd@0: transprob1 = 'ergodic';
wolffd@0: transprob2 = 'leftright';
wolffd@0: transprob3 = 'leftright';
wolffd@0: termprob2 = 'rightstop';
wolffd@0: termprob3 = 'rightstop';
wolffd@0: 
wolffd@0: 
wolffd@0: for i=1:2:nargs
wolffd@0:   switch args{i},
wolffd@0:    case 'discrete_obs', discrete_obs = args{i+1}; 
wolffd@0:    case 'Oargs',        Oargs = args{i+1};
wolffd@0:    case 'Q1args',       Q1args = args{i+1};
wolffd@0:    case 'Q2args',       Q2args = args{i+1};
wolffd@0:    case 'Q3args',       Q3args = args{i+1};
wolffd@0:    case 'F2args',       F2args = args{i+1};
wolffd@0:    case 'F3args',       F3args = args{i+1};
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: ns = zeros(1,ss);
wolffd@0: ns(Qnodes) = Qsizes;
wolffd@0: ns(obs) = Osize;
wolffd@0: ns(Fnodes) = 2;
wolffd@0: 
wolffd@0: dnodes = [Qnodes Fnodes];
wolffd@0: if discrete_obs
wolffd@0:   dnodes = [dnodes obs];
wolffd@0: end
wolffd@0: onodes = [obs];
wolffd@0: 
wolffd@0: bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'discrete', dnodes, 'names', names);
wolffd@0: eclass = bnet.equiv_class;
wolffd@0: 
wolffd@0: if strcmp(startprob1, 'ergodic')
wolffd@0:   startprob1 = normalise(ones(1,ns(Q1)));
wolffd@0: end
wolffd@0: if strcmp(startprob2, 'leftstart')
wolffd@0:   startprob2 = zeros(ns(Q1), ns(Q2));
wolffd@0:   starpbrob2(:, 1) = 1.0;
wolffd@0: end
wolffd@0: if strcmp(startprob3, 'leftstart')
wolffd@0:   startprob3 = zeros(ns(Q2), ns(Q3));
wolffd@0:   starpbrob3(:, 1) = 1.0;
wolffd@0: end
wolffd@0: 
wolffd@0: if strcmp(termprob2, 'rightstop')
wolffd@0:   p = 0.9;
wolffd@0:   termprob2 = zeros(Qsize(2),2);
wolffd@0:   termprob2(:, 2) = p; 
wolffd@0:   termprob2(:, 1) = 1-p; 
wolffd@0:   termprob2(1:(Qsize(2)-1), 1) = 1; 
wolffd@0: end
wolffd@0: if strcmp(termprob3, 'rightstop')
wolffd@0:   p = 0.9;
wolffd@0:   termprob3 = zeros(Qsize(3),2);
wolffd@0:   termprob3(:, 2) = p; 
wolffd@0:   termprob3(:, 1) = 1-p; 
wolffd@0:   termprob3(1:(Qsize(3)-1), 1) = 1; 
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % SLICE 1
wolffd@0: 
wolffd@0: % We clamp untied nodes in the first slice, since their params can't be estimated
wolffd@0: % from just one sequence
wolffd@0: 
wolffd@0: bnet.CPD{eclass(Q1,1)} = tabular_CPD(bnet, Q1, 'CPT', startprob1, 'adjustable', 0);
wolffd@0: bnet.CPD{eclass(Q2,1)} = tabular_CPD(bnet, Q2, 'CPT', startprob2, 'adjustable', 0);
wolffd@0: bnet.CPD{eclass(Q3,1)} = tabular_CPD(bnet, Q3, 'CPT', startprob3, 'adjustable', 0);
wolffd@0: 
wolffd@0: bnet.CPD{eclass(F2,1)}  = hhmmF_CPD(bnet, F2, Qnodes, 2, D, 'termprob', termprob2);
wolffd@0: bnet.CPD{eclass(F3,1)}  = tabular_CPD(bnet, F3, 'CPT', termprob3);
wolffd@0: 
wolffd@0: if discrete_obs
wolffd@0:   bnet.CPD{eclass(obs,1)} = tabular_CPD(bnet, obs, Oargs{:});
wolffd@0: else
wolffd@0:   bnet.CPD{eclass(obs,1)} = gaussian_CPD(bnet, obs, Oargs{:});
wolffd@0: end
wolffd@0: 
wolffd@0: % SLICE 2
wolffd@0: 
wolffd@0: bnet.CPD{eclass(Q1,2)} = hhmmQ_CPD(bnet, Q1+ss, Qnodes, 1, D, 'transprob', transprob1, 'startprob', startprob1);
wolffd@0: bnet.CPD{eclass(Q2,2)} = hhmmQ_CPD(bnet, Q2+ss, Qnodes, 2, D, 'transprob', transprob2, 'startprob', startprob2);
wolffd@0: bnet.CPD{eclass(Q3,2)} = hhmmQ_CPD(bnet, Q3+ss, Qnodes, 3, D, 'transprob', transprob3, 'startprob', startprob3);
wolffd@0: