diff options
Diffstat (limited to 'sourcecodes/bnt-master/HMM')
31 files changed, 1318 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/HMM/CVS/Entries b/sourcecodes/bnt-master/HMM/CVS/Entries new file mode 100644 index 00000000..b9105b50 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/CVS/Entries @@ -0,0 +1,29 @@ +/README.txt/1.1.1.1/Thu Jun 9 01:22:48 2005// +/dhmm_em.m/1.1.1.1/Thu Jun 9 01:25:04 2005// +/dhmm_em_demo.m/1.1.1.1/Sun May 4 22:01:12 2003// +/dhmm_em_online.m/1.1.1.1/Sun May 4 22:02:58 2003// +/dhmm_em_online_demo.m/1.1.1.1/Sun May 4 22:04:10 2003// +/dhmm_logprob.m/1.1.1.1/Sun May 4 22:01:34 2003// +/dhmm_logprob_brute_force.m/1.1.1.1/Wed May 29 15:59:56 2002// +/dhmm_logprob_path.m/1.1.1.1/Wed May 29 15:59:56 2002// +/dhmm_sample.m/1.1.1.1/Mon May 31 22:19:50 2004// +/dhmm_sample_endstate.m/1.1.1.1/Sun May 4 22:00:34 2003// +/fixed_lag_smoother.m/1.1.1.1/Wed Jan 22 17:56:04 2003// +/fixed_lag_smoother_demo.m/1.1.1.1/Thu Jun 9 01:27:20 2005// +/fwdback.m/1.1.1.1/Thu Jun 9 01:17:50 2005// +/gausshmm_train_observed.m/1.1.1.1/Thu Feb 12 23:08:22 2004// +/mc_sample.m/1.1.1.1/Mon May 24 22:26:34 2004// +/mc_sample_endstate.m/1.1.1.1/Wed Jan 22 20:32:28 2003// +/mdp_sample.m/1.1.1.1/Wed May 29 15:59:56 2002// +/mhmmParzen_train_observed.m/1.1.1.1/Sat Feb 14 02:06:30 2004// +/mhmm_em.m/1.1.1.1/Sun Feb 8 04:52:42 2004// +/mhmm_em_demo.m/1.1.1.1/Tue May 13 16:11:22 2003// +/mhmm_logprob.m/1.1.1.1/Sun May 4 22:11:54 2003// +/mhmm_sample.m/1.1.1.1/Wed May 26 00:32:28 2004// +/mk_leftright_transmat.m/1.1.1.1/Wed May 29 15:59:58 2002// +/mk_rightleft_transmat.m/1.1.1.1/Fri Nov 22 21:45:52 2002// +/pomdp_sample.m/1.1.1.1/Sun May 4 21:58:20 2003// +/testHMM.m/1.1.1.1/Thu Jun 9 01:25:50 2005// +/transmat_train_observed.m/1.1.1.1/Sun Aug 29 12:41:52 2004// +/viterbi_path.m/1.1.1.1/Sat Oct 23 01:18:22 2004// +D diff --git a/sourcecodes/bnt-master/HMM/CVS/Repository b/sourcecodes/bnt-master/HMM/CVS/Repository new file mode 100644 index 00000000..5f13f5c3 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/CVS/Repository @@ -0,0 +1 @@ +FullBNT/HMM diff --git a/sourcecodes/bnt-master/HMM/CVS/Root b/sourcecodes/bnt-master/HMM/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/HMM/README.txt b/sourcecodes/bnt-master/HMM/README.txt new file mode 100644 index 00000000..2e6ab78f --- /dev/null +++ b/sourcecodes/bnt-master/HMM/README.txt @@ -0,0 +1,23 @@ +Hidden Markov Model (HMM) Toolbox written by Kevin Murphy (1998). +See http://www.ai.mit.edu/~murphyk/Software/hmm.html for details. + +Models +------ + +dhmm = HMM with discrete output +mhmm = HMM with mixture of Gaussians output; + Use mhmm with M=1 components to simulate an HMM with a single Gaussian output. + +Demos +----- + +mhmm_em_demo +dhmm_em_demo +dhmm_em_online_demo +fixed_lag_smoother_demo + +References +----------- + +See "A tutorial on Hidden Markov Models and selected applications in speech recognition", + L. Rabiner, 1989, Proc. IEEE 77(2):257--286. diff --git a/sourcecodes/bnt-master/HMM/dhmm_em.m b/sourcecodes/bnt-master/HMM/dhmm_em.m new file mode 100644 index 00000000..34ac1c6d --- /dev/null +++ b/sourcecodes/bnt-master/HMM/dhmm_em.m @@ -0,0 +1,124 @@ +function [LL, prior, transmat, obsmat, nrIterations] = ... + dhmm_em(data, prior, transmat, obsmat, varargin) +% LEARN_DHMM Find the ML/MAP parameters of an HMM with discrete outputs using EM. +% [ll_trace, prior, transmat, obsmat, iterNr] = learn_dhmm(data, prior0, transmat0, obsmat0, ...) +% +% Notation: Q(t) = hidden state, Y(t) = observation +% +% INPUTS: +% data{ex} or data(ex,:) if all sequences have the same length +% prior(i) +% transmat(i,j) +% obsmat(i,o) +% +% Optional parameters may be passed as 'param_name', param_value pairs. +% Parameter names are shown below; default values in [] - if none, argument is mandatory. +% +% 'max_iter' - max number of EM iterations [10] +% 'thresh' - convergence threshold [1e-4] +% 'verbose' - if 1, print out loglik at every iteration [1] +% 'obs_prior_weight' - weight to apply to uniform dirichlet prior on observation matrix [0] +% +% To clamp some of the parameters, so learning does not change them: +% 'adj_prior' - if 0, do not change prior [1] +% 'adj_trans' - if 0, do not change transmat [1] +% 'adj_obs' - if 0, do not change obsmat [1] +% +% Modified by Herbert Jaeger so xi are not computed individually +% but only their sum (over time) as xi_summed; this is the only way how they are used +% and it saves a lot of memory. + +[max_iter, thresh, verbose, obs_prior_weight, adj_prior, adj_trans, adj_obs] = ... + process_options(varargin, 'max_iter', 10, 'thresh', 1e-4, 'verbose', 1, ... + 'obs_prior_weight', 0, 'adj_prior', 1, 'adj_trans', 1, 'adj_obs', 1); + +previous_loglik = -inf; +loglik = 0; +converged = 0; +num_iter = 1; +LL = []; + +if ~iscell(data) + data = num2cell(data, 2); % each row gets its own cell +end + +while (num_iter <= max_iter) & ~converged + % E step + [loglik, exp_num_trans, exp_num_visits1, exp_num_emit] = ... + compute_ess_dhmm(prior, transmat, obsmat, data, obs_prior_weight); + + % M step + if adj_prior + prior = normalise(exp_num_visits1); + end + if adj_trans & ~isempty(exp_num_trans) + transmat = mk_stochastic(exp_num_trans); + end + if adj_obs + obsmat = mk_stochastic(exp_num_emit); + end + + if verbose, fprintf(1, 'iteration %d, loglik = %f\n', num_iter, loglik); end + num_iter = num_iter + 1; + converged = em_converged(loglik, previous_loglik, thresh); + previous_loglik = loglik; + LL = [LL loglik]; +end +nrIterations = num_iter - 1; + +%%%%%%%%%%%%%%%%%%%%%%% + +function [loglik, exp_num_trans, exp_num_visits1, exp_num_emit, exp_num_visitsT] = ... + compute_ess_dhmm(startprob, transmat, obsmat, data, dirichlet) +% COMPUTE_ESS_DHMM Compute the Expected Sufficient Statistics for an HMM with discrete outputs +% function [loglik, exp_num_trans, exp_num_visits1, exp_num_emit, exp_num_visitsT] = ... +% compute_ess_dhmm(startprob, transmat, obsmat, data, dirichlet) +% +% INPUTS: +% startprob(i) +% transmat(i,j) +% obsmat(i,o) +% data{seq}(t) +% dirichlet - weighting term for uniform dirichlet prior on expected emissions +% +% OUTPUTS: +% exp_num_trans(i,j) = sum_l sum_{t=2}^T Pr(X(t-1) = i, X(t) = j| Obs(l)) +% exp_num_visits1(i) = sum_l Pr(X(1)=i | Obs(l)) +% exp_num_visitsT(i) = sum_l Pr(X(T)=i | Obs(l)) +% exp_num_emit(i,o) = sum_l sum_{t=1}^T Pr(X(t) = i, O(t)=o| Obs(l)) +% where Obs(l) = O_1 .. O_T for sequence l. + +numex = length(data); +[S O] = size(obsmat); +exp_num_trans = zeros(S,S); +exp_num_visits1 = zeros(S,1); +exp_num_visitsT = zeros(S,1); +exp_num_emit = dirichlet*ones(S,O); +loglik = 0; + +for ex=1:numex + obs = data{ex}; + T = length(obs); + %obslik = eval_pdf_cond_multinomial(obs, obsmat); + obslik = multinomial_prob(obs, obsmat); + [alpha, beta, gamma, current_ll, xi_summed] = fwdback(startprob, transmat, obslik); + + loglik = loglik + current_ll; + exp_num_trans = exp_num_trans + xi_summed; + exp_num_visits1 = exp_num_visits1 + gamma(:,1); + exp_num_visitsT = exp_num_visitsT + gamma(:,T); + % loop over whichever is shorter + if T < O + for t=1:T + o = obs(t); + exp_num_emit(:,o) = exp_num_emit(:,o) + gamma(:,t); + end + else + for o=1:O + ndx = find(obs==o); + if ~isempty(ndx) + exp_num_emit(:,o) = exp_num_emit(:,o) + sum(gamma(:, ndx), 2); + end + end + end +end diff --git a/sourcecodes/bnt-master/HMM/dhmm_em_demo.m b/sourcecodes/bnt-master/HMM/dhmm_em_demo.m new file mode 100644 index 00000000..c6ffa1dd --- /dev/null +++ b/sourcecodes/bnt-master/HMM/dhmm_em_demo.m @@ -0,0 +1,25 @@ +O = 3; +Q = 2; + +% "true" parameters +prior0 = normalise(rand(Q,1)); +transmat0 = mk_stochastic(rand(Q,Q)); +obsmat0 = mk_stochastic(rand(Q,O)); + +% training data +T = 5; +nex = 10; +data = dhmm_sample(prior0, transmat0, obsmat0, T, nex); + +% initial guess of parameters +prior1 = normalise(rand(Q,1)); +transmat1 = mk_stochastic(rand(Q,Q)); +obsmat1 = mk_stochastic(rand(Q,O)); + +% improve guess of parameters using EM +[LL, prior2, transmat2, obsmat2] = dhmm_em(data, prior1, transmat1, obsmat1, 'max_iter', 5); +LL + +% use model to compute log likelihood +loglik = dhmm_logprob(data, prior2, transmat2, obsmat2) +% log lik is slightly different than LL(end), since it is computed after the final M step diff --git a/sourcecodes/bnt-master/HMM/dhmm_em_online.m b/sourcecodes/bnt-master/HMM/dhmm_em_online.m new file mode 100644 index 00000000..605a1f5b --- /dev/null +++ b/sourcecodes/bnt-master/HMM/dhmm_em_online.m @@ -0,0 +1,80 @@ +function [transmat, obsmat, exp_num_trans, exp_num_emit, gamma, ll] = dhmm_em_online(... + prior, transmat, obsmat, exp_num_trans, exp_num_emit, decay, data, ... + act, adj_trans, adj_obs, dirichlet, filter_only) +% ONLINE_EM Adjust the parameters using a weighted combination of the old and new expected statistics +% +% [transmat, obsmat, exp_num_trans, exp_num_emit, gamma, ll] = online_em(... +% prior, transmat, obsmat, exp_num_trans, exp_num_emit, decay, data, act, ... +% adj_trans, adj_obs, dirichlet, filter_only) +% +% 0 < decay < 1, with smaller values meaning the past is forgotten more quickly. +% (We need to decay the old ess, since they were based on out-of-date parameters.) +% The other params are as in learn_hmm. +% We do a single forwards-backwards pass on the provided data, initializing with the specified prior. +% (If filter_only = 1, we only do a forwards pass.) + +if ~exist('act'), act = []; end +if ~exist('adj_trans'), adj_trans = 1; end +if ~exist('adj_obs'), adj_obs = 1; end +if ~exist('dirichlet'), dirichlet = 0; end +if ~exist('filter_only'), filter_only = 0; end + +% E step +olikseq = multinomial_prob(data, obsmat); +if isempty(act) + [alpha, beta, gamma, ll, xi] = fwdback(prior, transmat, olikseq, 'fwd_only', filter_only); +else + [alpha, beta, gamma, ll, xi] = fwdback(prior, transmat, olikseq, 'fwd_only', filter_only, ... + 'act', act); +end + +% Increment ESS +[S O] = size(obsmat); +if adj_obs + exp_num_emit = decay*exp_num_emit + dirichlet*ones(S,O); + T = length(data); + if T < O + for t=1:T + o = data(t); + exp_num_emit(:,o) = exp_num_emit(:,o) + gamma(:,t); + end + else + for o=1:O + ndx = find(data==o); + if ~isempty(ndx) + exp_num_emit(:,o) = exp_num_emit(:,o) + sum(gamma(:, ndx), 2); + end + end + end +end + +if adj_trans & (T > 1) + if isempty(act) + exp_num_trans = decay*exp_num_trans + sum(xi,3); + else + % act(2) determines Q(2), xi(:,:,1) holds P(Q(1), Q(2)) + A = length(transmat); + for a=1:A + ndx = find(act(2:end)==a); + if ~isempty(ndx) + exp_num_trans{a} = decay*exp_num_trans{a} + sum(xi(:,:,ndx), 3); + end + end + end +end + + +% M step + +if adj_obs + obsmat = mk_stochastic(exp_num_emit); +end +if adj_trans & (T>1) + if isempty(act) + transmat = mk_stochastic(exp_num_trans); + else + for a=1:A + transmat{a} = mk_stochastic(exp_num_trans{a}); + end + end +end diff --git a/sourcecodes/bnt-master/HMM/dhmm_em_online_demo.m b/sourcecodes/bnt-master/HMM/dhmm_em_online_demo.m new file mode 100644 index 00000000..579d89d5 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/dhmm_em_online_demo.m @@ -0,0 +1,93 @@ +% Example of online EM applied to a simple POMDP with fixed action seq + +clear all + +% Create a really easy model to learn +rand('state', 1); +O = 2; +S = 2; +A = 2; +prior0 = [1 0]'; +transmat0 = cell(1,A); +transmat0{1} = [0.9 0.1; 0.1 0.9]; % long runs of 1s and 2s +transmat0{2} = [0.1 0.9; 0.9 0.1]; % short runs +obsmat0 = eye(2); + +%prior0 = normalise(rand(S,1)); +%transmat0 = mk_stochastic(rand(S,S)); +%obsmat0 = mk_stochastic(rand(S,O)); + +T = 10; +act = [1*ones(1,25) 2*ones(1,25) 1*ones(1,25) 2*ones(1,25)]; +data = pomdp_sample(prior0, transmat0, obsmat0, act); +%data = sample_dhmm(prior0, transmat0, obsmat0, T, 1); + +% Initial guess of params +rand('state', 2); % different seed! +transmat1 = cell(1,A); +for a=1:A + transmat1{a} = mk_stochastic(rand(S,S)); +end +obsmat1 = mk_stochastic(rand(S,O)); +prior1 = prior0; % so it labels states the same way + +% Uniformative Dirichlet prior (expected sufficient statistics / pseudo counts) +e = 0.001; +ess_trans = cell(1,A); +for a=1:A + ess_trans{a} = repmat(e, S, S); +end +ess_emit = repmat(e, S, O); + +% Params +w = 2; +decay_sched = [0.1:0.1:0.9]; + +% Initialize +LL1 = zeros(1,T); +t = 1; +y = data(t); +data_win = y; +act_win = [1]; % arbitrary initial value +[prior1, LL1(1)] = normalise(prior1 .* obsmat1(:,y)); + +% Iterate +for t=2:T + y = data(t); + a = act(t); + if t <= w + data_win = [data_win y]; + act_win = [act_win a]; + else + data_win = [data_win(2:end) y]; + act_win = [act_win(2:end) a]; + prior1 = gamma(:, 2); + end + d = decay_sched(min(t, length(decay_sched))); + [transmat1, obsmat1, ess_trans, ess_emit, gamma, ll] = dhmm_em_online(... + prior1, transmat1, obsmat1, ess_trans, ess_emit, d, data_win, act_win); + bel = gamma(:, end); + LL1(t) = ll/length(data_win); + %fprintf('t=%d, ll=%f\n', t, ll); +end + +LL1(1) = LL1(2); % since initial likelihood is for 1 slice +plot(1:T, LL1, 'rx-'); + + +% compare with offline learning + +if 0 +rand('state', 2); % same seed as online learner +transmat2 = cell(1,A); +for a=1:A + transmat2{a} = mk_stochastic(rand(S,S)); +end +obsmat2 = mk_stochastic(rand(S,O)); +prior2 = prior0; +[LL2, prior2, transmat2, obsmat2] = dhmm_em(data, prior2, transmat2, obsmat2, .... + 'max_iter', 10, 'thresh', 1e-3, 'verbose', 1, 'act', act); + +LL2 = LL2 / T + +end diff --git a/sourcecodes/bnt-master/HMM/dhmm_logprob.m b/sourcecodes/bnt-master/HMM/dhmm_logprob.m new file mode 100644 index 00000000..f8afe835 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/dhmm_logprob.m @@ -0,0 +1,22 @@ +function [loglik, errors] = dhmm_logprob(data, prior, transmat, obsmat) +% LOG_LIK_DHMM Compute the log-likelihood of a dataset using a discrete HMM +% [loglik, errors] = log_lik_dhmm(data, prior, transmat, obsmat) +% +% data{m} or data(m,:) is the m'th sequence +% errors is a list of the cases which received a loglik of -infinity + +if ~iscell(data) + data = num2cell(data, 2); +end +ncases = length(data); + +loglik = 0; +errors = []; +for m=1:ncases + obslik = multinomial_prob(data{m}, obsmat); + [alpha, beta, gamma, ll] = fwdback(prior, transmat, obslik, 'fwd_only', 1); + if ll==-inf + errors = [errors m]; + end + loglik = loglik + ll; +end diff --git a/sourcecodes/bnt-master/HMM/dhmm_logprob_brute_force.m b/sourcecodes/bnt-master/HMM/dhmm_logprob_brute_force.m new file mode 100644 index 00000000..abffd0c9 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/dhmm_logprob_brute_force.m @@ -0,0 +1,21 @@ +function logp = enumerate_HMM_loglik(prior, transmat, obsmat) +% ENUMERATE_HMM_LOGLIK Compute the log likelihood of a sequence by exhaustive (O(Q^T)) enumeration. +% logp = enumerate_HMM_loglik(prior, transmat, obsmat) +% +% Inputs: +% prior(i) = Pr(Q(1) = i) +% transmat(i,j) = Pr(Q(t+1)=j | Q(t)=i) +% obsmat(i,t) = Pr(y(t) | Q(t)=i) + +Q = length(prior); +T = size(obsmat, 2); +sizes = repmat(Q, 1, T); + +psum = 0; +for i=1:Q^T + qs = ind2subv(sizes, i); % make the state sequence + psum = psum + prob_path(prior, transmat, obsmat, qs); +end +logp = log(psum) + + diff --git a/sourcecodes/bnt-master/HMM/dhmm_logprob_path.m b/sourcecodes/bnt-master/HMM/dhmm_logprob_path.m new file mode 100644 index 00000000..11a1e3af --- /dev/null +++ b/sourcecodes/bnt-master/HMM/dhmm_logprob_path.m @@ -0,0 +1,15 @@ +function [ll, p] = prob_path(prior, transmat, obsmat, qs) +% PROB_PATH Compute the prob. of a specific path (state sequence) through an HMM. +% [ll, p] = prob_path(prior, transmat, obsmat, states) +% +% ll = log prob path +% p(t) = Pr(O(t)) * Pr(Q(t) -> Q(t+1)) for 1<=t<T, p(T) = Pr(O(T)) + +T = size(obsmat, 2); +p = zeros(1,T); +p(1) = prior(qs(1)) * obsmat(qs(1),1); +for t=2:T + p(t) = transmat(qs(t-1), qs(t)) * obsmat(qs(t),t); +end + +ll = sum(log(p)); diff --git a/sourcecodes/bnt-master/HMM/dhmm_sample.m b/sourcecodes/bnt-master/HMM/dhmm_sample.m new file mode 100644 index 00000000..5122c81b --- /dev/null +++ b/sourcecodes/bnt-master/HMM/dhmm_sample.m @@ -0,0 +1,10 @@ +function [obs, hidden] = dhmm_sample(initial_prob, transmat, obsmat, numex, len) +% SAMPLE_DHMM Generate random sequences from a Hidden Markov Model with discrete outputs. +% +% [obs, hidden] = sample_dhmm(initial_prob, transmat, obsmat, numex, len) +% Each row of obs is an observation sequence of length len. + +hidden = mc_sample(initial_prob, transmat, len, numex); +obs = multinomial_sample(hidden, obsmat); + + diff --git a/sourcecodes/bnt-master/HMM/dhmm_sample_endstate.m b/sourcecodes/bnt-master/HMM/dhmm_sample_endstate.m new file mode 100644 index 00000000..1972ebd2 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/dhmm_sample_endstate.m @@ -0,0 +1,19 @@ +function [obs, hidden] = dhmm_sample_endstate(startprob, transmat, obsmat, endprob, numex) +% SAMPLE_DHMM Generate random sequences from an HMM with discrete outputs. +% function [obs, hidden] = sample_dhmm_endstate(startprob, transmat, obsmat, endprob, numex) +% +% We sample until we have have entered the end state +% obs{m} and hidden{m} are the m'th sequence + +hidden = cell(1,numex); +obs = cell(1,numex); + +for m=1:numex + hidden{m} = mc_sample_endstate(startprob, transmat, endprob); + T = length(hidden{m}); + obs{m} = zeros(1,T); + for t=1:T + h = hidden{m}(t); + obs{m}(t) = sample_discrete(obsmat(h,:)); + end +end diff --git a/sourcecodes/bnt-master/HMM/fixed_lag_smoother.m b/sourcecodes/bnt-master/HMM/fixed_lag_smoother.m new file mode 100644 index 00000000..9e19df98 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/fixed_lag_smoother.m @@ -0,0 +1,65 @@ +function [alpha, obslik, gamma, xi] = fixed_lag_smoother(d, alpha, obslik, obsvec, transmat, act) +% FIXED_LAG_SMOOTHER Computed smoothed posterior estimates within a window given previous filtered window. +% [alpha, obslik, gamma, xi] = fixed_lag_smoother(d, alpha, obslik, obsvec, transmat, act) +% +% d >= 2 is the desired window width. +% Actually, we use d=min(d, t0), where t0 is the current time. +% +% alpha(:, t0-d:t0-1) - length d window, excluding t0 (Columns indexed 1..d) +% obslik(:, t0-d:t0-1) - length d window +% obsvec - likelihood vector for current observation +% transmat - transition matrix +% If we specify the optional 'act' argument, transmat{a} should be a cell array, and +% act(t0-d:t0) - length d window, last column = current action +% +% Output: +% alpha(:, t0-d+1:t0) - last column = new filtered estimate +% obslik(:, t0-d+1:t0) - last column = obsvec +% xi(:, :, t0-d+1:t0-1) - 2 slice smoothed window +% gamma(:, t0-d+1:t0) - smoothed window +% +% As usual, we define (using T=d) +% alpha(i,t) = Pr(Q(t)=i | Y(1:t)) +% gamma(i,t) = Pr(Q(t)=i | Y(1:T)) +% xi(i,j,t) = Pr(Q(t)=i, Q(t+1)=j | Y(1:T)) +% +% obslik(i,t) = Pr(Y(t) | Q(t)=i) +% transmat{a}(i,j) = Pr(Q(t)=j | Q(t-1)=i, A(t)=a) + +[S n] = size(alpha); +d = min(d, n+1); +if d < 2 + error('must keep a window of length at least 2'); +end + +if ~exist('act') + act = ones(1, n+1); + transmat = { transmat }; +end + +% pluck out last d-1 components from the history +alpha = alpha(:, n-d+2:n); +obslik = obslik(:, n-d+2:n); + +% Extend window by 1 +t = d; +obslik(:,t) = obsvec; +xi = zeros(S, S, d-1); +xi(:,:,t-1) = normalise((alpha(:,t-1) * obslik(:,t)') .* transmat{act(t)}); +alpha(:,t) = sum(xi(:,:,t-1), 1)'; + +% Now smooth backwards inside the window +beta = ones(S, d); +T = d; +%fprintf('smooth from %d to 1, i.e., %d to %d\n', d, t0, t0-d+1); +gamma(:,T) = alpha(:,T); +for t=T-1:-1:1 + b = beta(:,t+1) .* obslik(:,t+1); + beta(:,t) = normalise(transmat{act(t)} * b); + gamma(:,t) = normalise(alpha(:,t) .* beta(:,t)); + xi(:,:,t) = normalise((transmat{act(t)} .* (alpha(:,t) * b'))); +end + + + + diff --git a/sourcecodes/bnt-master/HMM/fixed_lag_smoother_demo.m b/sourcecodes/bnt-master/HMM/fixed_lag_smoother_demo.m new file mode 100644 index 00000000..9cf19f90 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/fixed_lag_smoother_demo.m @@ -0,0 +1,38 @@ +% Example of fixed lag smoothing + +rand('state', 1); +S = 2; +O = 2; +T = 7; +data = sample_discrete([0.5 0.5], 1, T); +transmat = mk_stochastic(rand(S,S)); +obsmat = mk_stochastic(rand(S,O)); +obslik = multinomial_prob(data, obsmat); +prior = [0.5 0.5]'; + + +[alpha0, beta0, gamma0, ll0, xi0] = fwdback(prior, transmat, obslik); + +w = 3; +alpha1 = zeros(S, T); +gamma1 = zeros(S, T); +xi1 = zeros(S, S, T-1); +t = 1; +b = obsmat(:, data(t)); +olik_win = b; % window of conditional observation likelihoods +alpha_win = normalise(prior .* b); +alpha1(:,t) = alpha_win; +for t=2:T + [alpha_win, olik_win, gamma_win, xi_win] = ... + fixed_lag_smoother(w, alpha_win, olik_win, obsmat(:, data(t)), transmat); + alpha1(:,max(1,t-w+1):t) = alpha_win; + gamma1(:,max(1,t-w+1):t) = gamma_win; + xi1(:,:,max(1,t-w+1):t-1) = xi_win; +end + +e = 1e-1; +%assert(approxeq(alpha0, alpha1, e)); +assert(approxeq(gamma0(:, T-w+1:end), gamma1(:, T-w+1:end), e)); +%assert(approxeq(xi0(:,:,T-w+1:end), xi1(:,:,T-w+1:end), e)); + + diff --git a/sourcecodes/bnt-master/HMM/fwdback.m b/sourcecodes/bnt-master/HMM/fwdback.m new file mode 100644 index 00000000..d74d5490 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/fwdback.m @@ -0,0 +1,197 @@ +function [alpha, beta, gamma, loglik, xi_summed, gamma2] = fwdback(init_state_distrib, ... + transmat, obslik, varargin) +% FWDBACK Compute the posterior probs. in an HMM using the forwards backwards algo. +% +% [alpha, beta, gamma, loglik, xi, gamma2] = fwdback(init_state_distrib, transmat, obslik, ...) +% +% Notation: +% Y(t) = observation, Q(t) = hidden state, M(t) = mixture variable (for MOG outputs) +% A(t) = discrete input (action) (for POMDP models) +% +% INPUT: +% init_state_distrib(i) = Pr(Q(1) = i) +% transmat(i,j) = Pr(Q(t) = j | Q(t-1)=i) +% or transmat{a}(i,j) = Pr(Q(t) = j | Q(t-1)=i, A(t-1)=a) if there are discrete inputs +% obslik(i,t) = Pr(Y(t)| Q(t)=i) +% (Compute obslik using eval_pdf_xxx on your data sequence first.) +% +% Optional parameters may be passed as 'param_name', param_value pairs. +% Parameter names are shown below; default values in [] - if none, argument is mandatory. +% +% For HMMs with MOG outputs: if you want to compute gamma2, you must specify +% 'obslik2' - obslik(i,j,t) = Pr(Y(t)| Q(t)=i,M(t)=j) [] +% 'mixmat' - mixmat(i,j) = Pr(M(t) = j | Q(t)=i) [] +% +% For HMMs with discrete inputs: +% 'act' - act(t) = action performed at step t +% +% Optional arguments: +% 'fwd_only' - if 1, only do a forwards pass and set beta=[], gamma2=[] [0] +% 'scaled' - if 1, normalize alphas and betas to prevent underflow [1] +% 'maximize' - if 1, use max-product instead of sum-product [0] +% +% OUTPUTS: +% alpha(i,t) = p(Q(t)=i | y(1:t)) (or p(Q(t)=i, y(1:t)) if scaled=0) +% beta(i,t) = p(y(t+1:T) | Q(t)=i)*p(y(t+1:T)|y(1:t)) (or p(y(t+1:T) | Q(t)=i) if scaled=0) +% gamma(i,t) = p(Q(t)=i | y(1:T)) +% loglik = log p(y(1:T)) +% xi(i,j,t-1) = p(Q(t-1)=i, Q(t)=j | y(1:T)) - NO LONGER COMPUTED +% xi_summed(i,j) = sum_{t=}^{T-1} xi(i,j,t) - changed made by Herbert Jaeger +% gamma2(j,k,t) = p(Q(t)=j, M(t)=k | y(1:T)) (only for MOG outputs) +% +% If fwd_only = 1, these become +% alpha(i,t) = p(Q(t)=i | y(1:t)) +% beta = [] +% gamma(i,t) = p(Q(t)=i | y(1:t)) +% xi(i,j,t-1) = p(Q(t-1)=i, Q(t)=j | y(1:t)) +% gamma2 = [] +% +% Note: we only compute xi if it is requested as a return argument, since it can be very large. +% Similarly, we only compute gamma2 on request (and if using MOG outputs). +% +% Examples: +% +% [alpha, beta, gamma, loglik] = fwdback(pi, A, multinomial_prob(sequence, B)); +% +% [B, B2] = mixgauss_prob(data, mu, Sigma, mixmat); +% [alpha, beta, gamma, loglik, xi, gamma2] = fwdback(pi, A, B, 'obslik2', B2, 'mixmat', mixmat); + +if nargout >= 5, compute_xi = 1; else compute_xi = 0; end +if nargout >= 6, compute_gamma2 = 1; else compute_gamma2 = 0; end + +[obslik2, mixmat, fwd_only, scaled, act, maximize, compute_xi, compute_gamma2] = ... + process_options(varargin, ... + 'obslik2', [], 'mixmat', [], ... + 'fwd_only', 0, 'scaled', 1, 'act', [], 'maximize', 0, ... + 'compute_xi', compute_xi, 'compute_gamma2', compute_gamma2); + +[Q T] = size(obslik); + +if isempty(obslik2) + compute_gamma2 = 0; +end + +if isempty(act) + act = ones(1,T); + transmat = { transmat } ; +end + +scale = ones(1,T); + +% scale(t) = Pr(O(t) | O(1:t-1)) = 1/c(t) as defined by Rabiner (1989). +% Hence prod_t scale(t) = Pr(O(1)) Pr(O(2)|O(1)) Pr(O(3) | O(1:2)) ... = Pr(O(1), ... ,O(T)) +% or log P = sum_t log scale(t). +% Rabiner suggests multiplying beta(t) by scale(t), but we can instead +% normalise beta(t) - the constants will cancel when we compute gamma. + +loglik = 0; + +alpha = zeros(Q,T); +gamma = zeros(Q,T); +if compute_xi + xi_summed = zeros(Q,Q); +else + xi_summed = []; +end + +%%%%%%%%% Forwards %%%%%%%%%% + +t = 1; +alpha(:,1) = init_state_distrib(:) .* obslik(:,t); +if scaled + %[alpha(:,t), scale(t)] = normaliseC(alpha(:,t)); + [alpha(:,t), scale(t)] = normalise(alpha(:,t)); +end +assert(approxeq(sum(alpha(:,t)),1)) +for t=2:T + %trans = transmat(:,:,act(t-1))'; + trans = transmat{act(t-1)}; + if maximize + m = max_mult(trans', alpha(:,t-1)); + %A = repmat(alpha(:,t-1), [1 Q]); + %m = max(trans .* A, [], 1); + else + m = trans' * alpha(:,t-1); + end + alpha(:,t) = m(:) .* obslik(:,t); + if scaled + %[alpha(:,t), scale(t)] = normaliseC(alpha(:,t)); + [alpha(:,t), scale(t)] = normalise(alpha(:,t)); + end + if compute_xi & fwd_only % useful for online EM + %xi(:,:,t-1) = normaliseC((alpha(:,t-1) * obslik(:,t)') .* trans); + xi_summed = xi_summed + normalise((alpha(:,t-1) * obslik(:,t)') .* trans); + end + assert(approxeq(sum(alpha(:,t)),1)) +end +if scaled + if any(scale==0) + loglik = -inf; + else + loglik = sum(log(scale)); + end +else + loglik = log(sum(alpha(:,T))); +end + +if fwd_only + gamma = alpha; + beta = []; + gamma2 = []; + return; +end + +%%%%%%%%% Backwards %%%%%%%%%% + +beta = zeros(Q,T); +if compute_gamma2 + M = size(mixmat, 2); + gamma2 = zeros(Q,M,T); +else + gamma2 = []; +end + +beta(:,T) = ones(Q,1); +%gamma(:,T) = normaliseC(alpha(:,T) .* beta(:,T)); +gamma(:,T) = normalise(alpha(:,T) .* beta(:,T)); +t=T; +if compute_gamma2 + denom = obslik(:,t) + (obslik(:,t)==0); % replace 0s with 1s before dividing + gamma2(:,:,t) = obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M]) ./ repmat(denom, [1 M]); + %gamma2(:,:,t) = normaliseC(obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M])); % wrong! +end +for t=T-1:-1:1 + b = beta(:,t+1) .* obslik(:,t+1); + %trans = transmat(:,:,act(t)); + trans = transmat{act(t)}; + if maximize + B = repmat(b(:)', Q, 1); + beta(:,t) = max(trans .* B, [], 2); + else + beta(:,t) = trans * b; + end + if scaled + %beta(:,t) = normaliseC(beta(:,t)); + beta(:,t) = normalise(beta(:,t)); + end + %gamma(:,t) = normaliseC(alpha(:,t) .* beta(:,t)); + gamma(:,t) = normalise(alpha(:,t) .* beta(:,t)); + if compute_xi + %xi(:,:,t) = normaliseC((trans .* (alpha(:,t) * b'))); + xi_summed = xi_summed + normalise((trans .* (alpha(:,t) * b'))); + end + if compute_gamma2 + denom = obslik(:,t) + (obslik(:,t)==0); % replace 0s with 1s before dividing + gamma2(:,:,t) = obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M]) ./ repmat(denom, [1 M]); + %gamma2(:,:,t) = normaliseC(obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M])); + end +end + +% We now explain the equation for gamma2 +% Let zt=y(1:t-1,t+1:T) be all observations except y(t) +% gamma2(Q,M,t) = P(Qt,Mt|yt,zt) = P(yt|Qt,Mt,zt) P(Qt,Mt|zt) / P(yt|zt) +% = P(yt|Qt,Mt) P(Mt|Qt) P(Qt|zt) / P(yt|zt) +% Now gamma(Q,t) = P(Qt|yt,zt) = P(yt|Qt) P(Qt|zt) / P(yt|zt) +% hence +% P(Qt,Mt|yt,zt) = P(yt|Qt,Mt) P(Mt|Qt) [P(Qt|yt,zt) P(yt|zt) / P(yt|Qt)] / P(yt|zt) +% = P(yt|Qt,Mt) P(Mt|Qt) P(Qt|yt,zt) / P(yt|Qt) diff --git a/sourcecodes/bnt-master/HMM/gausshmm_train_observed.m b/sourcecodes/bnt-master/HMM/gausshmm_train_observed.m new file mode 100644 index 00000000..0e96b365 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/gausshmm_train_observed.m @@ -0,0 +1,44 @@ +function [initState, transmat, mu, Sigma] = gausshmm_train_observed(obsData, hiddenData, ... + nstates, varargin) +% GAUSSHMM_TRAIN_OBSERVED Estimate params of HMM with Gaussian output from fully observed sequences +% [initState, transmat, mu, Sigma] = gausshmm_train_observed(obsData, hiddenData, nstates,...) +% +% INPUT +% If all sequences have the same length +% obsData(:,t,ex) +% hiddenData(ex,t) - must be ROW vector if only one sequence +% If sequences have different lengths, we use cell arrays +% obsData{ex}(:,t) +% hiddenData{ex}(t) +% +% Optional argumnets +% dirichletPriorWeight - for smoothing transition matrix counts +% +% Optional parameters from mixgauss_Mstep: +% 'cov_type' - 'full', 'diag' or 'spherical' ['full'] +% 'tied_cov' - 1 (Sigma) or 0 (Sigma_i) [0] +% 'clamped_cov' - pass in clamped value, or [] if unclamped [ [] ] +% 'clamped_mean' - pass in clamped value, or [] if unclamped [ [] ] +% 'cov_prior' - Lambda_i, added to YY(:,:,i) [0.01*eye(d,d,Q)] +% +% Output +% mu(:,q) +% Sigma(:,:,q) + +[dirichletPriorWeight, other] = process_options(... + varargin, 'dirichletPriorWeight', 0); + +[transmat, initState] = transmat_train_observed(hiddenData, nstates, ... + 'dirichletPriorWeight', dirichletPriorWeight); + +% convert to obsData(:,t*nex) +if ~iscell(obsData) + [D T Nex] = size(obsData); + obsData = reshape(obsData, D, T*Nex); +else + obsData = cat(2, obsData{:}); + hiddenData = cat(2,hiddenData{:}); +end +[mu, Sigma] = condgaussTrainObserved(obsData, hiddenData(:), nstates, varargin{:}); + + diff --git a/sourcecodes/bnt-master/HMM/mc_sample.m b/sourcecodes/bnt-master/HMM/mc_sample.m new file mode 100644 index 00000000..87488f44 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/mc_sample.m @@ -0,0 +1,17 @@ +function S = mc_sample(prior, trans, len, numex) +% SAMPLE_MC Generate random sequences from a Markov chain. +% STATE = SAMPLE_MC(PRIOR, TRANS, LEN) generates a sequence of length LEN. +% +% STATE = SAMPLE_MC(PRIOR, TRANS, LEN, N) generates N rows each of length LEN. + +if nargin==3 + numex = 1; +end + +S = zeros(numex,len); +for i=1:numex + S(i, 1) = sample_discrete(prior); + for t=2:len + S(i, t) = sample_discrete(trans(S(i,t-1),:)); + end +end diff --git a/sourcecodes/bnt-master/HMM/mc_sample_endstate.m b/sourcecodes/bnt-master/HMM/mc_sample_endstate.m new file mode 100644 index 00000000..8dc9e0ba --- /dev/null +++ b/sourcecodes/bnt-master/HMM/mc_sample_endstate.m @@ -0,0 +1,27 @@ +function S = sample_mc_endstate(startprob, trans, endprob) +% SAMPLE_MC_ENDSTATE Generate a random sequence from a Markov chain until enter the endstate. +% seq = sample_mc(startprob, trans, endprob) + +% add an end state +Q = size(trans,1); +transprob = zeros(Q,Q+1); +end_state = Q+1; +for i=1:Q + for j=1:Q + transprob(i,j) = (1-endprob(i)) * trans(i,j); + end + transprob(i,end_state) = endprob(i); + %assert(approxeq(sum(transprob(i,:)), 1)) +end + +S = []; +S(1) = sample_discrete(startprob); +t = 1; +p = endprob(S(t)); +stop = (S(1) == end_state); +while ~stop + S(t+1) = sample_discrete(transprob(S(t),:)); + stop = (S(t+1) == end_state); + t = t + 1; +end +S = S(1:end-1); % don't include end state diff --git a/sourcecodes/bnt-master/HMM/mdp_sample.m b/sourcecodes/bnt-master/HMM/mdp_sample.m new file mode 100644 index 00000000..63923c9c --- /dev/null +++ b/sourcecodes/bnt-master/HMM/mdp_sample.m @@ -0,0 +1,18 @@ +function state = sample_mdp(prior, trans, act) +% SAMPLE_MDP Sample a sequence of states from a Markov Decision Process. +% state = sample_mdp(prior, trans, act) +% +% Inputs: +% prior(i) = Pr(Q(1)=i) +% trans{a}(i,j) = Pr(Q(t)=j | Q(t-1)=i, A(t)=a) +% act(a) = A(t), so act(1) is ignored +% +% Output: +% state is a vector of length T=length(act) + +len = length(act); +state = zeros(1,len); +state(1) = sample_discrete(prior); +for t=2:len + state(t) = sample_discrete(trans{act(t)}(state(t-1),:)); +end diff --git a/sourcecodes/bnt-master/HMM/mhmmParzen_train_observed.m b/sourcecodes/bnt-master/HMM/mhmmParzen_train_observed.m new file mode 100644 index 00000000..1786e9a2 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/mhmmParzen_train_observed.m @@ -0,0 +1,37 @@ +function [initState, transmat, mu, Nproto, pick] = mhmmParzen_train_observed(obsData, hiddenData, ... + nstates, maxNproto, varargin) +% mhmmParzentrain_observed with mixture of Gaussian outputs from fully observed sequences +% function [initState, transmat, mu, Nproto] = mhmm_train_observed_parzen(obsData, hiddenData, ... +% nstates, maxNproto) +% +% +% INPUT +% If all sequences have the same length +% obsData(:,t,ex) +% hiddenData(ex,t) - must be ROW vector if only one sequence +% If sequences have different lengths, we use cell arrays +% obsData{ex}(:,t) +% hiddenData{ex}(t) +% +% Optional argumnets +% dirichletPriorWeight - for smoothing transition matrix counts +% mkSymmetric +% +% Output +% mu(:,q) +% Nproto(q) is the number of prototypes (mixture components) chosen for state q + +[transmat, initState] = transmat_train_observed(... + hiddenData, nstates, varargin{:}); + +% convert to obsData(:,t*nex) +if ~iscell(obsData) + [D T Nex] = size(obsData); + obsData = reshape(obsData, D, T*Nex); +else + obsData = cat(2, obsData{:}); + hiddenData = cat(2, hiddenData{:}); +end +[mu, Nproto, pick] = parzen_fit_select_unif(obsData, hiddenData(:), maxNproto); + + diff --git a/sourcecodes/bnt-master/HMM/mhmm_em.m b/sourcecodes/bnt-master/HMM/mhmm_em.m new file mode 100644 index 00000000..036c5222 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/mhmm_em.m @@ -0,0 +1,173 @@ +function [LL, prior, transmat, mu, Sigma, mixmat] = ... + mhmm_em(data, prior, transmat, mu, Sigma, mixmat, varargin); +% LEARN_MHMM Compute the ML parameters of an HMM with (mixtures of) Gaussians output using EM. +% [ll_trace, prior, transmat, mu, sigma, mixmat] = learn_mhmm(data, ... +% prior0, transmat0, mu0, sigma0, mixmat0, ...) +% +% Notation: Q(t) = hidden state, Y(t) = observation, M(t) = mixture variable +% +% INPUTS: +% data{ex}(:,t) or data(:,t,ex) if all sequences have the same length +% prior(i) = Pr(Q(1) = i), +% transmat(i,j) = Pr(Q(t+1)=j | Q(t)=i) +% mu(:,j,k) = E[Y(t) | Q(t)=j, M(t)=k ] +% Sigma(:,:,j,k) = Cov[Y(t) | Q(t)=j, M(t)=k] +% mixmat(j,k) = Pr(M(t)=k | Q(t)=j) : set to [] or ones(Q,1) if only one mixture component +% +% Optional parameters may be passed as 'param_name', param_value pairs. +% Parameter names are shown below; default values in [] - if none, argument is mandatory. +% +% 'max_iter' - max number of EM iterations [10] +% 'thresh' - convergence threshold [1e-4] +% 'verbose' - if 1, print out loglik at every iteration [1] +% 'cov_type' - 'full', 'diag' or 'spherical' ['full'] +% +% To clamp some of the parameters, so learning does not change them: +% 'adj_prior' - if 0, do not change prior [1] +% 'adj_trans' - if 0, do not change transmat [1] +% 'adj_mix' - if 0, do not change mixmat [1] +% 'adj_mu' - if 0, do not change mu [1] +% 'adj_Sigma' - if 0, do not change Sigma [1] +% +% If the number of mixture components differs depending on Q, just set the trailing +% entries of mixmat to 0, e.g., 2 components if Q=1, 3 components if Q=2, +% then set mixmat(1,3)=0. In this case, B2(1,3,:)=1.0. + +if ~isstr(varargin{1}) % catch old syntax + error('optional arguments should be passed as string/value pairs') +end + +[max_iter, thresh, verbose, cov_type, adj_prior, adj_trans, adj_mix, adj_mu, adj_Sigma] = ... + process_options(varargin, 'max_iter', 10, 'thresh', 1e-4, 'verbose', 1, ... + 'cov_type', 'full', 'adj_prior', 1, 'adj_trans', 1, 'adj_mix', 1, ... + 'adj_mu', 1, 'adj_Sigma', 1); + +previous_loglik = -inf; +loglik = 0; +converged = 0; +num_iter = 1; +LL = []; + +if ~iscell(data) + data = num2cell(data, [1 2]); % each elt of the 3rd dim gets its own cell +end +numex = length(data); + + +O = size(data{1},1); +Q = length(prior); +if isempty(mixmat) + mixmat = ones(Q,1); +end +M = size(mixmat,2); +if M == 1 + adj_mix = 0; +end + +while (num_iter <= max_iter) & ~converged + % E step + [loglik, exp_num_trans, exp_num_visits1, postmix, m, ip, op] = ... + ess_mhmm(prior, transmat, mixmat, mu, Sigma, data); + + + % M step + if adj_prior + prior = normalise(exp_num_visits1); + end + if adj_trans + transmat = mk_stochastic(exp_num_trans); + end + if adj_mix + mixmat = mk_stochastic(postmix); + end + if adj_mu | adj_Sigma + [mu2, Sigma2] = mixgauss_Mstep(postmix, m, op, ip, 'cov_type', cov_type); + if adj_mu + mu = reshape(mu2, [O Q M]); + end + if adj_Sigma + Sigma = reshape(Sigma2, [O O Q M]); + end + end + + if verbose, fprintf(1, 'iteration %d, loglik = %f\n', num_iter, loglik); end + num_iter = num_iter + 1; + converged = em_converged(loglik, previous_loglik, thresh); + previous_loglik = loglik; + LL = [LL loglik]; +end + + +%%%%%%%%% + +function [loglik, exp_num_trans, exp_num_visits1, postmix, m, ip, op] = ... + ess_mhmm(prior, transmat, mixmat, mu, Sigma, data) +% ESS_MHMM Compute the Expected Sufficient Statistics for a MOG Hidden Markov Model. +% +% Outputs: +% exp_num_trans(i,j) = sum_l sum_{t=2}^T Pr(Q(t-1) = i, Q(t) = j| Obs(l)) +% exp_num_visits1(i) = sum_l Pr(Q(1)=i | Obs(l)) +% +% Let w(i,k,t,l) = P(Q(t)=i, M(t)=k | Obs(l)) +% where Obs(l) = Obs(:,:,l) = O_1 .. O_T for sequence l +% Then +% postmix(i,k) = sum_l sum_t w(i,k,t,l) (posterior mixing weights/ responsibilities) +% m(:,i,k) = sum_l sum_t w(i,k,t,l) * Obs(:,t,l) +% ip(i,k) = sum_l sum_t w(i,k,t,l) * Obs(:,t,l)' * Obs(:,t,l) +% op(:,:,i,k) = sum_l sum_t w(i,k,t,l) * Obs(:,t,l) * Obs(:,t,l)' + + +verbose = 0; + +%[O T numex] = size(data); +numex = length(data); +O = size(data{1},1); +Q = length(prior); +M = size(mixmat,2); +exp_num_trans = zeros(Q,Q); +exp_num_visits1 = zeros(Q,1); +postmix = zeros(Q,M); +m = zeros(O,Q,M); +op = zeros(O,O,Q,M); +ip = zeros(Q,M); + +mix = (M>1); + +loglik = 0; +if verbose, fprintf(1, 'forwards-backwards example # '); end +for ex=1:numex + if verbose, fprintf(1, '%d ', ex); end + %obs = data(:,:,ex); + obs = data{ex}; + T = size(obs,2); + if mix + [B, B2] = mixgauss_prob(obs, mu, Sigma, mixmat); + [alpha, beta, gamma, current_loglik, xi, gamma2] = ... + fwdback(prior, transmat, B, 'obslik2', B2, 'mixmat', mixmat); + else + B = mixgauss_prob(obs, mu, Sigma); + [alpha, beta, gamma, current_loglik, xi] = fwdback(prior, transmat, B); + end + loglik = loglik + current_loglik; + if verbose, fprintf(1, 'll at ex %d = %f\n', ex, loglik); end + + exp_num_trans = exp_num_trans + sum(xi,3); + exp_num_visits1 = exp_num_visits1 + gamma(:,1); + + if mix + postmix = postmix + sum(gamma2,3); + else + postmix = postmix + sum(gamma,2); + gamma2 = reshape(gamma, [Q 1 T]); % gamma2(i,m,t) = gamma(i,t) + end + for i=1:Q + for k=1:M + w = reshape(gamma2(i,k,:), [1 T]); % w(t) = w(i,k,t,l) + wobs = obs .* repmat(w, [O 1]); % wobs(:,t) = w(t) * obs(:,t) + m(:,i,k) = m(:,i,k) + sum(wobs, 2); % m(:) = sum_t w(t) obs(:,t) + op(:,:,i,k) = op(:,:,i,k) + wobs * obs'; % op(:,:) = sum_t w(t) * obs(:,t) * obs(:,t)' + ip(i,k) = ip(i,k) + sum(sum(wobs .* obs, 2)); % ip = sum_t w(t) * obs(:,t)' * obs(:,t) + end + end +end +if verbose, fprintf(1, '\n'); end diff --git a/sourcecodes/bnt-master/HMM/mhmm_em_demo.m b/sourcecodes/bnt-master/HMM/mhmm_em_demo.m new file mode 100644 index 00000000..3e361545 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/mhmm_em_demo.m @@ -0,0 +1,40 @@ +if 1 + O = 4; + T = 10; + nex = 50; + M = 2; + Q = 3; +else + O = 8; %Number of coefficients in a vector + T = 420; %Number of vectors in a sequence + nex = 1; %Number of sequences + M = 1; %Number of mixtures + Q = 6; %Number of states +end +cov_type = 'full'; + +data = randn(O,T,nex); + +% initial guess of parameters +prior0 = normalise(rand(Q,1)); +transmat0 = mk_stochastic(rand(Q,Q)); + +if 0 + Sigma0 = repmat(eye(O), [1 1 Q M]); + % Initialize each mean to a random data point + indices = randperm(T*nex); + mu0 = reshape(data(:,indices(1:(Q*M))), [O Q M]); + mixmat0 = mk_stochastic(rand(Q,M)); +else + [mu0, Sigma0] = mixgauss_init(Q*M, data, cov_type); + mu0 = reshape(mu0, [O Q M]); + Sigma0 = reshape(Sigma0, [O O Q M]); + mixmat0 = mk_stochastic(rand(Q,M)); +end + +[LL, prior1, transmat1, mu1, Sigma1, mixmat1] = ... + mhmm_em(data, prior0, transmat0, mu0, Sigma0, mixmat0, 'max_iter', 5); + + +loglik = mhmm_logprob(data, prior1, transmat1, mu1, Sigma1, mixmat1); + diff --git a/sourcecodes/bnt-master/HMM/mhmm_logprob.m b/sourcecodes/bnt-master/HMM/mhmm_logprob.m new file mode 100644 index 00000000..e27a4eb9 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/mhmm_logprob.m @@ -0,0 +1,30 @@ +function [loglik, errors] = mhmm_logprob(data, prior, transmat, mu, Sigma, mixmat) +% LOG_LIK_MHMM Compute the log-likelihood of a dataset using a (mixture of) Gaussians HMM +% [loglik, errors] = log_lik_mhmm(data, prior, transmat, mu, sigma, mixmat) +% +% data{m}(:,t) or data(:,t,m) if all cases have same length +% errors is a list of the cases which received a loglik of -infinity +% +% Set mixmat to ones(Q,1) or omit it if there is only 1 mixture component + +Q = length(prior); +if size(mixmat,1) ~= Q % trap old syntax + error('mixmat should be QxM') +end +if nargin < 6, mixmat = ones(Q,1); end + +if ~iscell(data) + data = num2cell(data, [1 2]); % each elt of the 3rd dim gets its own cell +end +ncases = length(data); + +loglik = 0; +errors = []; +for m=1:ncases + obslik = mixgauss_prob(data{m}, mu, Sigma, mixmat); + [alpha, beta, gamma, ll] = fwdback(prior, transmat, obslik, 'fwd_only', 1); + if ll==-inf + errors = [errors m]; + end + loglik = loglik + ll; +end diff --git a/sourcecodes/bnt-master/HMM/mhmm_sample.m b/sourcecodes/bnt-master/HMM/mhmm_sample.m new file mode 100644 index 00000000..83034cd4 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/mhmm_sample.m @@ -0,0 +1,31 @@ +function [obs, hidden] = mhmm_sample(T, numex, initial_prob, transmat, mu, Sigma, mixmat) +% SAMPLE_MHMM Generate random sequences from an HMM with (mixtures of) Gaussian output. +% [obs, hidden] = sample_mhmm(T, numex, initial_prob, transmat, mu, Sigma, mixmat) +% +% INPUTS: +% T - length of each sequence +% numex - num. sequences +% init_state_prob(i) = Pr(Q(1) = i) +% transmat(i,j) = Pr(Q(t+1)=j | Q(t)=i) +% mu(:,j,k) = mean of Y(t) given Q(t)=j, M(t)=k +% Sigma(:,:,j,k) = cov. of Y(t) given Q(t)=j, M(t)=k +% mixmat(j,k) = Pr(M(t)=k | Q(t)=j) : set to ones(Q,1) or omit if single mixture +% +% OUTPUT: +% obs(:,t,l) = observation vector at time t for sequence l +% hidden(t,l) = the hidden state at time t for sequence l + +Q = length(initial_prob); +if nargin < 7, mixmat = ones(Q,1); end +O = size(mu,1); +hidden = zeros(T, numex); +obs = zeros(O, T, numex); + +hidden = mc_sample(initial_prob, transmat, T, numex)'; +for i=1:numex + for t=1:T + q = hidden(t,i); + m = sample_discrete(mixmat(q,:), 1, 1); + obs(:,t,i) = gaussian_sample(mu(:,q,m), Sigma(:,:,q,m), 1); + end +end diff --git a/sourcecodes/bnt-master/HMM/mk_leftright_transmat.m b/sourcecodes/bnt-master/HMM/mk_leftright_transmat.m new file mode 100644 index 00000000..8a36e1a6 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/mk_leftright_transmat.m @@ -0,0 +1,6 @@ +function transmat = mk_leftright_transmat(Q, p) +% MK_LEFTRIGHT_TRANSMAT Q = num states, p = prob on (i,i), 1-p on (i,i+1) +% function transmat = mk_leftright_transmat(Q, p) + +transmat = p*diag(ones(Q,1)) + (1-p)*diag(ones(Q-1,1),1); +transmat(Q,Q)=1; diff --git a/sourcecodes/bnt-master/HMM/mk_rightleft_transmat.m b/sourcecodes/bnt-master/HMM/mk_rightleft_transmat.m new file mode 100644 index 00000000..ce7cf62b --- /dev/null +++ b/sourcecodes/bnt-master/HMM/mk_rightleft_transmat.m @@ -0,0 +1,6 @@ +function transmat = mk_rightleft_transmat(Q, p) +% MK_RIGHTLEFT_TRANSMAT Q = num states, p = prob on (i,i), 1-p on (i,i+1) +% function transmat = mk_rightleft_transmat(Q, p) + +transmat = p*diag(ones(Q,1)) + (1-p)*diag(ones(Q-1,1),-1); +transmat(1,1)=1; diff --git a/sourcecodes/bnt-master/HMM/pomdp_sample.m b/sourcecodes/bnt-master/HMM/pomdp_sample.m new file mode 100644 index 00000000..836ce39f --- /dev/null +++ b/sourcecodes/bnt-master/HMM/pomdp_sample.m @@ -0,0 +1,20 @@ +function [obs, hidden] = pomdp_sample(initial_prob, transmat, obsmat, act) +% SAMPLE_POMDP Generate a random sequence from a Partially Observed Markov Decision Process. +% [obs, hidden] = sample_pomdp(prior, transmat, obsmat, act) +% +% Inputs: +% prior(i) = Pr(Q(1)=i) +% transmat{a}(i,j) = Pr(Q(t)=j | Q(t-1)=i, A(t)=a) +% obsmat(i,k) = Pr(Y(t)=k | Q(t)=i) +% act(a) = A(t), so act(1) is ignored +% +% Output: +% obs and hidden are vectors of length T=length(act) + + +len = length(act); +hidden = mdp_sample(initial_prob, transmat, act); +obs = zeros(1, len); +for t=1:len + obs(t) = sample_discrete(obsmat(hidden(t),:)); +end diff --git a/sourcecodes/bnt-master/HMM/testHMM.m b/sourcecodes/bnt-master/HMM/testHMM.m new file mode 100644 index 00000000..b3ced724 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/testHMM.m @@ -0,0 +1,5 @@ +% Run all the demos, to check everything is "syntactically correct" +mhmm_em_demo +dhmm_em_demo +dhmm_em_online_demo +fixed_lag_smoother_demo diff --git a/sourcecodes/bnt-master/HMM/transmat_train_observed.m b/sourcecodes/bnt-master/HMM/transmat_train_observed.m new file mode 100644 index 00000000..06918429 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/transmat_train_observed.m @@ -0,0 +1,39 @@ +function [transmat, initState] = transmat_train_observed(labels, nstates, varargin) +% transmat_train_observed ML estimation from fully observed data +% function [transmat, initState] = transmat_train_observed(labels, nstates, varargin) +% +% If all sequences have the same length +% labels(ex,t) +% If sequences have different lengths, we use cell arrays +% labels{ex}(t) + +[dirichletPriorWeight, mkSymmetric, other] = process_options(... + varargin, 'dirichletPriorWeight', 0, 'mkSymmetric', 0); + +if ~iscell(labels) + [numex T] = size(labels); + if T==1 + labels = labels'; + end + %fprintf('T=%d, numex=%d\n', T, numex); + labels = num2cell(labels,2); % each row gets its own cell +end +numex = length(labels); + +counts = zeros(nstates, nstates); +counts1 = zeros(nstates,1); +for s=1:numex + labs = labels{s}; labs = labs(:)'; + dat = [labs(1:end-1); labs(2:end)]; + counts = counts + compute_counts(dat, [nstates nstates]); + q = labs(1); + counts1(q) = counts1(q) + 1; +end +pseudo_counts = dirichletPriorWeight*ones(nstates, nstates); +if mkSymmetric + counts = counts + counts'; +end +transmat = mk_stochastic(counts + pseudo_counts); +initState = normalize(counts1 + dirichletPriorWeight*ones(nstates,1)); + + diff --git a/sourcecodes/bnt-master/HMM/viterbi_path.m b/sourcecodes/bnt-master/HMM/viterbi_path.m new file mode 100644 index 00000000..88010277 --- /dev/null +++ b/sourcecodes/bnt-master/HMM/viterbi_path.m @@ -0,0 +1,62 @@ +function path = viterbi_path(prior, transmat, obslik) +% VITERBI Find the most-probable (Viterbi) path through the HMM state trellis. +% path = viterbi(prior, transmat, obslik) +% +% Inputs: +% prior(i) = Pr(Q(1) = i) +% transmat(i,j) = Pr(Q(t+1)=j | Q(t)=i) +% obslik(i,t) = Pr(y(t) | Q(t)=i) +% +% Outputs: +% path(t) = q(t), where q1 ... qT is the argmax of the above expression. + + +% delta(j,t) = prob. of the best sequence of length t-1 and then going to state j, and O(1:t) +% psi(j,t) = the best predecessor state, given that we ended up in state j at t + +scaled = 1; + +T = size(obslik, 2); +prior = prior(:); +Q = length(prior); + +delta = zeros(Q,T); +psi = zeros(Q,T); +path = zeros(1,T); +scale = ones(1,T); + + +t=1; +delta(:,t) = prior .* obslik(:,t); +if scaled + [delta(:,t), n] = normalise(delta(:,t)); + scale(t) = 1/n; +end +psi(:,t) = 0; % arbitrary value, since there is no predecessor to t=1 +for t=2:T + for j=1:Q + [delta(j,t), psi(j,t)] = max(delta(:,t-1) .* transmat(:,j)); + delta(j,t) = delta(j,t) * obslik(j,t); + end + if scaled + [delta(:,t), n] = normalise(delta(:,t)); + scale(t) = 1/n; + end +end +[p, path(T)] = max(delta(:,T)); +for t=T-1:-1:1 + path(t) = psi(path(t+1),t+1); +end + +% If scaled==0, p = prob_path(best_path) +% If scaled==1, p = Pr(replace sum with max and proceed as in the scaled forwards algo) +% Both are different from p(data) as computed using the sum-product (forwards) algorithm + +if 0 +if scaled + loglik = -sum(log(scale)); + %loglik = prob_path(prior, transmat, obslik, path); +else + loglik = log(p); +end +end |