Mercurial > hg > camir-aes2014
comparison toolboxes/FullBNT-1.0.7/bnt/examples/dynamic/HHMM/Motif/learn_motif_hhmm.m @ 0:e9a9cd732c1e tip
first hg version after svn
author | wolffd |
---|---|
date | Tue, 10 Feb 2015 15:05:51 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:e9a9cd732c1e |
---|---|
1 | |
2 seed = 0; | |
3 rand('state', seed); | |
4 randn('state', seed); | |
5 | |
6 chars = ['a', 'c', 'g', 't']; | |
7 motif = 'accca'; | |
8 motif_length = length(motif); | |
9 motif_code = zeros(1, motif_length); | |
10 for i=1:motif_length | |
11 motif_code(i) = find(chars == motif(i)); | |
12 end | |
13 | |
14 [bnet_init, Qnodes, Fnodes, Onode] = mk_motif_hhmm('motif_length', length(motif)); | |
15 %[bnet_init, Qnodes, Fnodes, Onode] = mk_motif_hhmm('motif_pattern', motif); | |
16 ss = bnet_init.nnodes_per_slice; | |
17 | |
18 | |
19 | |
20 % We generate a training set by creating uniform sequences, | |
21 % and inserting a single motif at a random location. | |
22 ntrain = 100; | |
23 T = 20; | |
24 cases = cell(1, ntrain); | |
25 | |
26 if 1 | |
27 % uniform background | |
28 background_dist = normalise(ones(1, length(chars))); | |
29 end | |
30 if 0 | |
31 % use a constant background | |
32 background_dist = zeros(1, length(chars)); | |
33 m = find(chars=='t'); | |
34 background_dist(m) = 1.0; | |
35 end | |
36 if 0 | |
37 % use a background skewed away from the motif | |
38 p = 0.01; q = (1-(2*p))/2; | |
39 background_dist = [p p q q]; | |
40 end | |
41 | |
42 unif_pos = normalise(ones(1, T-length(motif))); | |
43 cases = cell(1, ntrain); | |
44 data = zeros(1,T); | |
45 for i=1:ntrain | |
46 data = sample_discrete(background_dist, 1, T); | |
47 L = sample_discrete(unif_pos, 1, 1); | |
48 data(L:L+length(motif)-1) = motif_code; | |
49 cases{i} = cell(ss, T); | |
50 cases{i}(Onode,:) = num2cell(data); | |
51 end | |
52 disp('sample training cases') | |
53 for i=1:5 | |
54 chars(cell2num(cases{i}(Onode,:))) | |
55 end | |
56 | |
57 engine_init = hmm_inf_engine(bnet_init); | |
58 | |
59 [bnet_learned, LL, engine_learned] = ... | |
60 learn_params_dbn_em(engine_init, cases, 'max_iter', 100, 'thresh', 1e-2); | |
61 % 'anneal', 1, 'anneal_rate', 0.7); | |
62 | |
63 % extract the learned motif profile | |
64 eclass = bnet_learned.equiv_class; | |
65 CPDO=struct(bnet_learned.CPD{eclass(Onode,1)}); | |
66 fprintf('columns = chars, rows = states\n'); | |
67 profile_learned = squeeze(CPDO.CPT(2,:,:)) | |
68 [m,ndx] = max(profile_learned, [], 2); | |
69 map_motif_learned = chars(ndx) | |
70 back_learned = squeeze(CPDO.CPT(1,1,:))' | |
71 %map_back_learned = chars(argmax(back_learned)) | |
72 | |
73 CPDO_init = struct(bnet_init.CPD{eclass(Onode,1)}); | |
74 profile_init = squeeze(CPDO_init.CPT(2,:,:)); | |
75 back_init = squeeze(CPDO_init.CPT(1,1,:))'; |