wolffd@0: function obslik = mk_hmm_obs_lik_vec(engine, evidence)
wolffd@0: 
wolffd@0: % P(o1,o2| h) = P(o1|h) * P(o2|h) where h = Q1,Q2,...
wolffd@0: 
wolffd@0: bnet = bnet_from_engine(engine);
wolffd@0: ss = length(bnet.intra);
wolffd@0: onodes = bnet.observed;
wolffd@0: hnodes = mysetdiff(1:ss, onodes);
wolffd@0: ns = bnet.node_sizes(:);
wolffd@0: ns(onodes) = 1;
wolffd@0: 
wolffd@0: Q = length(engine.startprob);
wolffd@0: obslik = ones(Q, 1);
wolffd@0: 
wolffd@0: for i=1:length(onodes)
wolffd@0:   o = onodes(i);
wolffd@0:   %data = cell2num(evidence(o,1));
wolffd@0:   data = evidence{o,1};
wolffd@0:   if myismember(o, bnet.dnodes)
wolffd@0:     obslik_i = eval_pdf_cond_multinomial(data, engine.obsprob{i}.CPT);
wolffd@0:   else
wolffd@0:     if bnet.auto_regressive(o)
wolffd@0:       error('can''t handle AR nodes')
wolffd@0:     end
wolffd@0:     %% calling mk_ghmm_obs_lik, which calls gaussian_prob, is slow, so we inline it
wolffd@0:     %% and use the pre-computed  inverse matrix
wolffd@0:     %obslik_i = mk_ghmm_obs_lik(data, engine.obsprob{i}.mu, engine.obsprob{i}.Sigma);
wolffd@0:     x = data(:);
wolffd@0:     m = engine.obsprob{i}.mu;
wolffd@0:     Qi = size(m, 2);
wolffd@0:     obslik_i = size(Qi, 1);
wolffd@0:     invC = engine.obsprob{i}.inv_Sigma;
wolffd@0:     denom = engine.obsprob{i}.denom;
wolffd@0:     for j=1:Qi
wolffd@0:       numer = exp(-0.5 * (x-m(:,j))' * invC(:,:,j) * (x-m(:,j)));
wolffd@0:       obslik_i(j) = numer / denom(j);
wolffd@0:     end
wolffd@0:   end
wolffd@0:   % convert P(o|ps) into P(o|h) by multiplying onto a (h,o) potential of all 1s
wolffd@0:   ps = bnet.parents{o};
wolffd@0:   dom = [ps o];
wolffd@0:   obspot_i = dpot(dom, ns(dom), obslik_i);
wolffd@0:   dom = [hnodes o];
wolffd@0:   obspot = dpot(dom, ns(dom));
wolffd@0:   obspot = multiply_by_pot(obspot, obspot_i);
wolffd@0:   % compute p(oi|h) * p(oj|h)
wolffd@0:   S = struct(obspot);
wolffd@0:   obslik = obslik .* S.T(:);
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: