wolffd@0: % Try to learn a 3 level HHMM similar to mk_square_hhmm
wolffd@0: % from synthetic discrete sequences
wolffd@0: 
wolffd@0: 
wolffd@0: discrete_obs = 1;
wolffd@0: supervised = 0;
wolffd@0: obs_finalF2 = 0;
wolffd@0: 
wolffd@0: seed = 1;
wolffd@0: rand('state', seed);
wolffd@0: randn('state', seed);
wolffd@0: 
wolffd@0: bnet_init = mk_square_hhmm(discrete_obs, 0);
wolffd@0: 
wolffd@0: ss = 6;
wolffd@0: Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;
wolffd@0: Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];
wolffd@0: 
wolffd@0: if supervised
wolffd@0:   bnet_init.observed = [Q1 Q2 Onode];
wolffd@0: else
wolffd@0:   bnet_init.observed = [Onode];
wolffd@0: end
wolffd@0: 
wolffd@0: if obs_finalF2
wolffd@0:   engine_init = jtree_dbn_inf_engine(bnet_init);
wolffd@0:   % can't use ndx version because sometimes F2 is hidden, sometimes observed
wolffd@0:   error('can''t observe F when learning')
wolffd@0:   % It is not possible to observe F2 if we learn
wolffd@0:   % because the update_ess method for hhmmF_CPD and hhmmQ_CPD assume
wolffd@0:   % the F nodes are always hidden (for speed).
wolffd@0:   % However, for generating, we might want to set the final F2=true
wolffd@0:   % to force all subroutines to finish.
wolffd@0: else
wolffd@0:   if supervised
wolffd@0:     engine_init = jtree_ndx_dbn_inf_engine(bnet_init);
wolffd@0:   else
wolffd@0:     engine_init = hmm_inf_engine(bnet_init);
wolffd@0:   end
wolffd@0: end
wolffd@0:   
wolffd@0: % generate some synthetic data (easier to debug)
wolffd@0: chars = ['L', 'l', 'U', 'u', 'R', 'r', 'D', 'd'];
wolffd@0: L=find(chars=='L'); l=find(chars=='l');
wolffd@0: U=find(chars=='U'); u=find(chars=='u');
wolffd@0: R=find(chars=='R'); r=find(chars=='r');
wolffd@0: D=find(chars=='D'); d=find(chars=='d');
wolffd@0: 
wolffd@0: cases = {};
wolffd@0: 
wolffd@0: T = 8;
wolffd@0: ev = cell(ss, T);
wolffd@0: ev(Onode,:) = num2cell([L l U u R r D d]);
wolffd@0: if supervised
wolffd@0:   ev(Q1,:) = num2cell(1*ones(1,T));
wolffd@0:   ev(Q2,:) = num2cell( [1 1 2 2 3 3 4 4]);
wolffd@0: end
wolffd@0: cases{1} = ev;
wolffd@0: cases{3} = ev;
wolffd@0: 
wolffd@0: T  = 8;
wolffd@0: ev = cell(ss, T);
wolffd@0: %we start with R then r, even though we are running the model 'backwards'!
wolffd@0: ev(Onode,:) = num2cell([R r U u L l D d]);
wolffd@0: 
wolffd@0: if supervised
wolffd@0:   ev(Q1,:) = num2cell(2*ones(1,T));
wolffd@0:   ev(Q2,:) = num2cell( [3 3 2 2 1 1 4 4]);
wolffd@0: end
wolffd@0: 
wolffd@0: cases{2} = ev;
wolffd@0: cases{4} = ev;
wolffd@0: 
wolffd@0: if obs_finalF2
wolffd@0:   for i=1:length(cases)
wolffd@0:     T = size(cases{i},2);
wolffd@0:     cases{i}(F2,T)={2}; % force F2 to be finished at end of seq
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % startprob should be shared for t=1:T,
wolffd@0: % but in the DBN it is shared for t=2:T,
wolffd@0: % so we train using a single long sequence.
wolffd@0: long_seq = cat(2, cases{:});
wolffd@0: [bnet_learned, LL, engine_learned] = ...
wolffd@0:     learn_params_dbn_em(engine_init, {long_seq}, 'max_iter', 200);
wolffd@0: 
wolffd@0: % figure out which subsequence each model is responsible for
wolffd@0: mpe = calc_mpe_dbn(engine_learned, long_seq);
wolffd@0: pretty_print_hhmm_parse(mpe, Qnodes, Fnodes, Onode, chars);
wolffd@0: 
wolffd@0: 
wolffd@0: % The "true" segmentation of the training sequence  is
wolffd@0: % Q1: 1                 2
wolffd@0: % O:  L l U u R r D d | R r U u L l D d | etc.
wolffd@0: % 
wolffd@0: % When we learn in a supervised fashion, we recover the "truth".
wolffd@0: 
wolffd@0: % When we learn in an unsupervised fashion with seed=1, we get
wolffd@0: % Q1: 2                       1
wolffd@0: % O:  L l U u R r D d  R r | U u L l D d | etc.
wolffd@0: %
wolffd@0: % This means for model 1:
wolffd@0: % starts in state 2
wolffd@0: % transitions 2->1, 1->4, 4->e, 3->2
wolffd@0: %
wolffd@0: % For model 2,
wolffd@0: % starts in state 1
wolffd@0: % transitions 1->2, 2->3, 3->4 or e, 4->3
wolffd@0: 
wolffd@0: % examine the params
wolffd@0: eclass = bnet_learned.equiv_class;
wolffd@0: CPDQ1=struct(bnet_learned.CPD{eclass(Q1,2)});
wolffd@0: CPDQ2=struct(bnet_learned.CPD{eclass(Q2,2)});
wolffd@0: CPDQ3=struct(bnet_learned.CPD{eclass(Q3,2)});
wolffd@0: CPDF2=struct(bnet_learned.CPD{eclass(F2,1)});
wolffd@0: CPDF3=struct(bnet_learned.CPD{eclass(F3,1)});
wolffd@0: CPDO=struct(bnet_learned.CPD{eclass(Onode,1)});
wolffd@0: 
wolffd@0: A_learned =add_hhmm_end_state(CPDQ2.transprob, CPDF2.termprob(:,:,2));
wolffd@0: squeeze(A_learned(:,1,:))
wolffd@0: squeeze(A_learned(:,2,:))
wolffd@0: 
wolffd@0: 
wolffd@0: % Does the "true" model have higher likelihood than the learned one?
wolffd@0: % i.e., Does the unsupervised method learn the wrong model because
wolffd@0: % we have the wrong cost fn, or because of local minima?
wolffd@0: 
wolffd@0: bnet_true = mk_square_hhmm(discrete_obs,1);
wolffd@0: 
wolffd@0: % examine the params
wolffd@0: eclass = bnet_learned.equiv_class;
wolffd@0: CPDQ1_true=struct(bnet_true.CPD{eclass(Q1,2)});
wolffd@0: CPDQ2_true=struct(bnet_true.CPD{eclass(Q2,2)});
wolffd@0: CPDQ3_true=struct(bnet_true.CPD{eclass(Q3,2)});
wolffd@0: CPDF2_true=struct(bnet_true.CPD{eclass(F2,1)});
wolffd@0: CPDF3_true=struct(bnet_true.CPD{eclass(F3,1)});
wolffd@0: 
wolffd@0: A_true =add_hhmm_end_state(CPDQ2_true.transprob, CPDF2_true.termprob(:,:,2));
wolffd@0: squeeze(A_true(:,1,:))
wolffd@0: 
wolffd@0: 
wolffd@0: if supervised
wolffd@0:   engine_true = jtree_ndx_dbn_inf_engine(bnet_true);
wolffd@0: else
wolffd@0:   engine_true = hmm_inf_engine(bnet_true);
wolffd@0: end
wolffd@0: 
wolffd@0: %[engine_learned, ll_learned] = enter_evidence(engine_learned, long_seq);
wolffd@0: %[engine_true, ll_true] = enter_evidence(engine_true, long_seq);
wolffd@0: [engine_learned, ll_learned] = enter_evidence(engine_learned, cases{2});
wolffd@0: [engine_true, ll_true] = enter_evidence(engine_true, cases{2});
wolffd@0: ll_learned
wolffd@0: ll_true
wolffd@0: 
wolffd@0: 
wolffd@0: % remove concatentation artefacts
wolffd@0: ll_learned = 0;
wolffd@0: ll_true = 0;
wolffd@0: for m=1:length(cases)
wolffd@0:   [engine_learned, ll_learned_tmp] = enter_evidence(engine_learned, cases{m});
wolffd@0:   [engine_true, ll_true_tmp] = enter_evidence(engine_true, cases{m});
wolffd@0:   ll_learned = ll_learned + ll_learned_tmp;
wolffd@0:   ll_true = ll_true + ll_true_tmp;
wolffd@0: end
wolffd@0: ll_learned
wolffd@0: ll_true
wolffd@0: 
wolffd@0: % In both cases, ll_learned >> ll_true
wolffd@0: % which shows we are using the wrong cost function!