package aima.core.probability.hmm.exact;
import java.util.LinkedList;
import java.util.List;
import aima.core.probability.CategoricalDistribution;
import aima.core.probability.hmm.HiddenMarkovModel;
import aima.core.probability.proposition.AssignmentProposition;
import aima.core.util.math.Matrix;
/**
* Artificial Intelligence A Modern Approach (3rd Edition): page 580.<br>
* <br>
*
* <pre>
* function FIXED-LAG-SMOOTHING(e<sub>t</sub>, hmm, d) returns a distribution over <b>X</b><sub>t-d</sub>
* inputs: e<sub>t</sub>, the current evidence from time step t
* hmm, a hidden Markov model with S * S transition matrix <b>T</b>
* d, the length of the lag for smoothing
* persistent: t, the current time, initially 1
* <b>f</b>, the forward message <b>P</b>(X<sub>t</sub> | e<sub>1:t</sub>), initially hmm.PRIOR
* <b>B</b>, the d-step backward transformation matrix, initially the identity matrix
* e<sub>t-d:t</sub>, double-ended list of evidence from t-d to t, initially empty
* local variables: <b>O</b><sub>t-d</sub>, <b>O</b><sub>t</sub>, diagonal matrices containing the sensor model information
*
* add e<sub>t</sub> to the end of e<sub>t-d:t</sub>
* <b>O</b><sub>t</sub> <- diagonal matrix containing <b>P</b>(e<sub>t</sub> | X<sub>t</sub>)
* if t > d then
* <b>f</b> <- FORWARD(<b>f</b>, e<sub>t</sub>)
* remove e<sub>t-d-1</sub> from the beginning of e<sub>t-d:t</sub>
* <b>O</b><sub>t-d</sub> <- diagonal matrix containing <b>P</b>(e<sub>t-d</sub> | X<sub>t-d</sub>)
* <b>B</b> <- <b>O</b><sup>-1</sup><sub>t-d</sub><b>B</b><b>T</b><b>O</b><sub>t</sub>
* else <b>B</b> <- <b>BTO</b><sub>t</sub>
* t <- t + 1
* if t > d then return NORMALIZE(<b>f</b> * <b>B1</b>) else return null
* </pre>
*
* Figure 15.6 An algorithm for smoothing with a fixed time lag of d steps,
* implemented as an online algorithm that outputs the new smoothed estimate
* given the observation for a new time step. Notice that the final output
* NORMALIZE(<b>f</b> * <b>B1</b>) is just α<b>f</b>*<b>b</b>, by Equation
* (15.14).<br>
* <br>
* <b>Note:</b> There appears to be two minor defects in the algorithm outlined
* in the book:<br>
* <b>f</b> <- FORWARD(<b>f</b>, e<sub>t</sub>)<br>
* should be:<br>
* <b>f</b> <- FORWARD(<b>f</b>, e<sub>t-d</sub>)<br>
* as we are returning a smoothed step for t-d and not the current time t. <br>
* <br>
* The update of:<br>
* t <- t + 1<br>
* should occur after the return value is calculated. Otherwise when t == d the
* value returned is based on HMM.prior in the calculation as opposed to a
* correctly calculated forward message. Comments welcome.
*
* @author Ciaran O'Reilly
* @author Ravi Mohan
*
*/
public class FixedLagSmoothing {
// persistent:
// t, the current time, initially 1
private int t = 1;
// <b>f</b>, the forward message <b>P</b>(X<sub>t</sub> | e<sub>1:t</sub>),
// initially hmm.PRIOR
private Matrix f = null;
// <b>B</b>, the d-step backward transformation matrix, initially the
// identity matrix
private Matrix B = null;
// e<sub>t-d:t</sub>, double-ended list of evidence from t-d to t, initially
// empty
private List<Matrix> e_tmd_to_t = new LinkedList<Matrix>();
// a hidden Markov model with S * S transition matrix <b>T</b>
private HiddenMarkovModel hmm = null;
// d, the length of the lag for smoothing
private int d = 1;
//
private Matrix unitMessage = null;
/**
* Create a Fixed-Lag-Smoothing implementation, that sets up the required
* persistent values.
*
* @param hmm
* a hidden Markov model with S * S transition matrix <b>T</b>
* @param d
* d, the length of the lag for smoothing
*/
public FixedLagSmoothing(HiddenMarkovModel hmm, int d) {
this.hmm = hmm;
this.d = d;
initPersistent();
}
/**
* Algorithm for smoothing with a fixed time lag of d steps, implemented as
* an online algorithm that outputs the new smoothed estimate given the
* observation for a new time step.
*
* @param et
* the current evidence from time step t
* @return a distribution over <b>X</b><sub>t-d</sub>
*/
public CategoricalDistribution fixedLagSmoothing(
List<AssignmentProposition> et) {
// local variables: <b>O</b><sub>t-d</sub>, <b>O</b><sub>t</sub>,
// diagonal matrices containing the sensor model information
Matrix O_tmd, O_t;
// add e<sub>t</sub> to the end of e<sub>t-d:t</sub>
e_tmd_to_t.add(hmm.getEvidence(et));
// <b>O</b><sub>t</sub> <- diagonal matrix containing
// <b>P</b>(e<sub>t</sub> | X<sub>t</sub>)
O_t = e_tmd_to_t.get(e_tmd_to_t.size() - 1);
// if t > d then
if (t > d) {
// remove e<sub>t-d-1</sub> from the beginning of e<sub>t-d:t</sub>
e_tmd_to_t.remove(0);
// <b>O</b><sub>t-d</sub> <- diagonal matrix containing
// <b>P</b>(e<sub>t-d</sub> | X<sub>t-d</sub>)
O_tmd = e_tmd_to_t.get(0);
// <b>f</b> <- FORWARD(<b>f</b>, e<sub>t-d</sub>)
f = forward(f, O_tmd);
// <b>B</b> <-
// <b>O</b><sup>-1</sup><sub>t-d</sub><b>B</b><b>T</b><b>O</b><sub>t</sub>
B = O_tmd.inverse().times(hmm.getTransitionModel().inverse())
.times(B).times(hmm.getTransitionModel()).times(O_t);
} else {
// else <b>B</b> <- <b>BTO</b><sub>t</sub>
B = B.times(hmm.getTransitionModel()).times(O_t);
}
// if t > d then return NORMALIZE(<b>f</b> * <b>B1</b>) else return null
CategoricalDistribution rVal = null;
if (t > d) {
rVal = hmm
.convert(hmm.normalize(f.arrayTimes(B.times(unitMessage))));
}
// t <- t + 1
t = t + 1;
return rVal;
}
/**
* The forward equation (15.5) in Matrix form becomes (15.12):<br>
*
* <pre>
* <b>f</b><sub>1:t+1</sub> = α<b>O</b><sub>t+1</sub><b>T</b><sup>T</sup><b>f</b><sub>1:t</sub>
* </pre>
*
* @param f1_t
* <b>f</b><sub>1:t</sub>
* @param O_tp1
* <b>O</b><sub>t+1</sub>
* @return <b>f</b><sub>1:t+1</sub>
*/
public Matrix forward(Matrix f1_t, Matrix O_tp1) {
return hmm.normalize(O_tp1.times(hmm.getTransitionModel().transpose()
.times(f1_t)));
}
//
// PRIVATE METHODS
//
private void initPersistent() {
// t, the current time, initially 1
t = 1;
// <b>f</b>, the forward message <b>P</b>(X<sub>t</sub> |
// e<sub>1:t</sub>),
// initially hmm.PRIOR
f = hmm.getPrior();
// <b>B</b>, the d-step backward transformation matrix, initially the
// identity matrix
B = Matrix.identity(f.getRowDimension(), f.getRowDimension());
// e<sub>t-d:t</sub>, double-ended list of evidence from t-d to t,
// initially
// empty
e_tmd_to_t.clear();
unitMessage = hmm.createUnitMessage();
}
}