|
wolffd@0
|
1 % Make the following HHMM
|
|
wolffd@0
|
2 %
|
|
wolffd@0
|
3 % S1 <----------------------> S2
|
|
wolffd@0
|
4 % | |
|
|
wolffd@0
|
5 % | |
|
|
wolffd@0
|
6 % M1 -> M2 -> M3 -> end B1 -> end
|
|
wolffd@0
|
7 %
|
|
wolffd@0
|
8 % where Mi represents the i'th letter in the motif
|
|
wolffd@0
|
9 % and B is the background state.
|
|
wolffd@0
|
10 % Si chooses between running the motif or the background.
|
|
wolffd@0
|
11 % The Si and B states have self loops (not shown).
|
|
wolffd@0
|
12
|
|
wolffd@0
|
13 if 0
|
|
wolffd@0
|
14 seed = 0;
|
|
wolffd@0
|
15 rand('state', seed);
|
|
wolffd@0
|
16 randn('state', seed);
|
|
wolffd@0
|
17 end
|
|
wolffd@0
|
18
|
|
wolffd@0
|
19 chars = ['a', 'c', 'g', 't'];
|
|
wolffd@0
|
20 Osize = length(chars);
|
|
wolffd@0
|
21
|
|
wolffd@0
|
22 motif_pattern = 'acca';
|
|
wolffd@0
|
23 motif_length = length(motif_pattern);
|
|
wolffd@0
|
24 Qsize = [2 motif_length];
|
|
wolffd@0
|
25 Qnodes = 1:2;
|
|
wolffd@0
|
26 D = 2;
|
|
wolffd@0
|
27 transprob = cell(1,D);
|
|
wolffd@0
|
28 termprob = cell(1,D);
|
|
wolffd@0
|
29 startprob = cell(1,D);
|
|
wolffd@0
|
30
|
|
wolffd@0
|
31 % startprob{d}(k,j), startprob{1}(1,j)
|
|
wolffd@0
|
32 % transprob{d}(i,k,j), transprob{1}(i,j)
|
|
wolffd@0
|
33 % termprob{d}(k,j)
|
|
wolffd@0
|
34
|
|
wolffd@0
|
35
|
|
wolffd@0
|
36 % LEVEL 1
|
|
wolffd@0
|
37
|
|
wolffd@0
|
38 startprob{1} = zeros(1, 2);
|
|
wolffd@0
|
39 startprob{1} = [1 0]; % always start in the background model
|
|
wolffd@0
|
40
|
|
wolffd@0
|
41 % When in the background state, we stay there with high prob
|
|
wolffd@0
|
42 % When in the motif state, we immediately return to the background state.
|
|
wolffd@0
|
43 transprob{1} = [0.8 0.2;
|
|
wolffd@0
|
44 1.0 0.0];
|
|
wolffd@0
|
45
|
|
wolffd@0
|
46
|
|
wolffd@0
|
47 % LEVEL 2
|
|
wolffd@0
|
48 startprob{2} = 'leftstart'; % both submodels start in substate 1
|
|
wolffd@0
|
49 transprob{2} = zeros(motif_length, 2, motif_length);
|
|
wolffd@0
|
50 termprob{2} = zeros(2, motif_length);
|
|
wolffd@0
|
51
|
|
wolffd@0
|
52 % In the background model, we only use state 1.
|
|
wolffd@0
|
53 transprob{2}(1,1,1) = 1; % self loop
|
|
wolffd@0
|
54 termprob{2}(1,1) = 0.2; % prob transition to end state
|
|
wolffd@0
|
55
|
|
wolffd@0
|
56 % Motif model
|
|
wolffd@0
|
57 transprob{2}(:,2,:) = mk_leftright_transmat(motif_length, 0);
|
|
wolffd@0
|
58 termprob{2}(2,end) = 1.0; % last state immediately terminates
|
|
wolffd@0
|
59
|
|
wolffd@0
|
60
|
|
wolffd@0
|
61 % OBS LEVEl
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63 obsprob = zeros([Qsize Osize]);
|
|
wolffd@0
|
64 if 0
|
|
wolffd@0
|
65 % uniform background model
|
|
wolffd@0
|
66 obsprob(1,1,:) = normalise(ones(Osize,1));
|
|
wolffd@0
|
67 else
|
|
wolffd@0
|
68 % deterministic background model (easy to see!)
|
|
wolffd@0
|
69 m = find(chars=='t');
|
|
wolffd@0
|
70 obsprob(1,1,m) = 1.0;
|
|
wolffd@0
|
71 end
|
|
wolffd@0
|
72 if 1
|
|
wolffd@0
|
73 % initialise with true motif (cheating)
|
|
wolffd@0
|
74 for i=1:motif_length
|
|
wolffd@0
|
75 m = find(chars == motif_pattern(i));
|
|
wolffd@0
|
76 obsprob(2,i,m) = 1.0;
|
|
wolffd@0
|
77 end
|
|
wolffd@0
|
78 end
|
|
wolffd@0
|
79
|
|
wolffd@0
|
80 Oargs = {'CPT', obsprob};
|
|
wolffd@0
|
81
|
|
wolffd@0
|
82 [bnet, Qnodes, Fnodes, Onode] = mk_hhmm('Qsizes', Qsize, 'Osize', Osize, 'discrete_obs', 1, ...
|
|
wolffd@0
|
83 'Oargs', Oargs, 'Ops', Qnodes(1:2), ...
|
|
wolffd@0
|
84 'startprob', startprob, 'transprob', transprob, 'termprob', termprob);
|
|
wolffd@0
|
85
|
|
wolffd@0
|
86
|
|
wolffd@0
|
87 Tmax = 20;
|
|
wolffd@0
|
88 usecell = 0;
|
|
wolffd@0
|
89
|
|
wolffd@0
|
90 for seqi=1:5
|
|
wolffd@0
|
91 evidence = sample_dbn(bnet, Tmax, usecell);
|
|
wolffd@0
|
92 chars(evidence(end,:))
|
|
wolffd@0
|
93 %T = size(evidence, 2)
|
|
wolffd@0
|
94 %pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, chars);
|
|
wolffd@0
|
95 end
|
