wolffd@0: % Make the following HHMM
wolffd@0: %
wolffd@0: %    S1 <----------------------> S2
wolffd@0: %    |                           |
wolffd@0: %    |                           |
wolffd@0: %   M1 -> M2 -> M3 -> end        B1 -> end
wolffd@0: %
wolffd@0: % where Mi represents the i'th letter in the motif
wolffd@0: % and B is the background state.
wolffd@0: % Si chooses between running the motif or the background.
wolffd@0: % The Si and B states have self loops (not shown).
wolffd@0: 
wolffd@0: if 0
wolffd@0: seed = 0;
wolffd@0: rand('state', seed);
wolffd@0: randn('state', seed);
wolffd@0: end
wolffd@0: 
wolffd@0: chars = ['a', 'c', 'g', 't'];
wolffd@0: Osize = length(chars);
wolffd@0: 
wolffd@0: motif_pattern = 'acca';
wolffd@0: motif_length = length(motif_pattern);
wolffd@0: Qsize = [2 motif_length];
wolffd@0: Qnodes = 1:2;
wolffd@0: D = 2;
wolffd@0: transprob = cell(1,D);
wolffd@0: termprob = cell(1,D);
wolffd@0: startprob = cell(1,D);
wolffd@0: 
wolffd@0: % startprob{d}(k,j), startprob{1}(1,j)
wolffd@0: % transprob{d}(i,k,j), transprob{1}(i,j)
wolffd@0: % termprob{d}(k,j)
wolffd@0: 
wolffd@0: 
wolffd@0: % LEVEL 1
wolffd@0: 
wolffd@0: startprob{1} = zeros(1, 2);
wolffd@0: startprob{1} = [1 0]; % always start in the background model
wolffd@0: 
wolffd@0: % When in the background state, we stay there with high prob
wolffd@0: % When in the motif state, we immediately return to the background state.
wolffd@0: transprob{1} = [0.8 0.2;
wolffd@0: 		1.0 0.0];
wolffd@0: 
wolffd@0: 
wolffd@0: % LEVEL 2
wolffd@0: startprob{2} = 'leftstart'; % both submodels start in substate 1
wolffd@0: transprob{2} = zeros(motif_length, 2, motif_length);
wolffd@0: termprob{2} = zeros(2, motif_length);
wolffd@0: 
wolffd@0: % In the background model, we only use state 1.
wolffd@0: transprob{2}(1,1,1) = 1; % self loop
wolffd@0: termprob{2}(1,1) = 0.2; % prob transition to end state
wolffd@0: 
wolffd@0: % Motif model
wolffd@0: transprob{2}(:,2,:) = mk_leftright_transmat(motif_length, 0);
wolffd@0: termprob{2}(2,end) = 1.0; % last state immediately terminates
wolffd@0: 
wolffd@0: 
wolffd@0: % OBS LEVEl
wolffd@0: 
wolffd@0: obsprob = zeros([Qsize Osize]);
wolffd@0: if 0
wolffd@0:   % uniform background model
wolffd@0:   obsprob(1,1,:) = normalise(ones(Osize,1));
wolffd@0: else
wolffd@0:   % deterministic background model (easy to see!)
wolffd@0:   m = find(chars=='t');
wolffd@0:   obsprob(1,1,m) = 1.0;
wolffd@0: end
wolffd@0: if 1
wolffd@0:   % initialise with true motif (cheating)
wolffd@0:   for i=1:motif_length
wolffd@0:     m = find(chars == motif_pattern(i));
wolffd@0:     obsprob(2,i,m) = 1.0;
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: Oargs = {'CPT', obsprob};
wolffd@0: 
wolffd@0: [bnet, Qnodes, Fnodes, Onode] = mk_hhmm('Qsizes', Qsize, 'Osize', Osize, 'discrete_obs', 1, ...
wolffd@0: 	       'Oargs', Oargs, 'Ops', Qnodes(1:2), ...
wolffd@0: 	       'startprob', startprob, 'transprob', transprob, 'termprob', termprob);
wolffd@0: 
wolffd@0: 
wolffd@0: Tmax = 20;
wolffd@0: usecell = 0;
wolffd@0: 
wolffd@0: for seqi=1:5
wolffd@0:   evidence = sample_dbn(bnet, Tmax, usecell);
wolffd@0:   chars(evidence(end,:))
wolffd@0:   %T = size(evidence, 2)
wolffd@0:   %pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, chars);
wolffd@0: end