Merge pull request #2142 from aabghari/master

HMM: calculate likelihood for data stream with/without pre-calculated emission probability
This commit is contained in:
Ryan Curtin
2021-01-12 08:58:25 -05:00
committed by GitHub
4 changed files with 505 additions and 31 deletions
+3
View File
@@ -12,6 +12,9 @@
* Add Triplet Margin Loss function (#2762).
* Add finalizers to Julia binding model types to fix memory handling (#2756).
* HMM: add functions to calculate likelihood for data stream with/without
pre-calculated emission probability (#2142).
### mlpack 3.4.2
###### 2020-10-26
+99 -2
View File
@@ -290,6 +290,79 @@ class HMM
*/
double LogLikelihood(const arma::mat& dataSeq) const;
/**
* Compute the log of the scaling factor of the given emission probability
* at time t. To calculate the log-likelihood for the whole sequence,
* accumulate log scale over the entire sequence
* This is meant for incremental or streaming computation of the
* log-likelihood of a sequence. For the first data point, provide an empty
* forwardLogProb vector.
*
* @param emissionLogProb emission probability at time t.
* @param forwardLogProb Vector in which forward probabilities will be saved.
* Passing forwardLogProb as an empty vector indicates the start of the
* sequence (i.e. time t=0).
* @return Log scale factor of the given sequence of emission at time t.
*/
double EmissionLogScaleFactor(const arma::vec& emissionLogProb,
arma::vec& forwardLogProb) const;
/**
* Compute the log-likelihood of the given emission probability up to time t,
* storing the result in logLikelihood.
* This is meant for incremental or streaming computation of the
* log-likelihood of a sequence. For the first data point, provide an empty
* forwardLogProb vector.
*
* @param emissionLogProb emission probability at time t.
* @param logLikelihood Log-likelihood of the given sequence of emission
* probability up to time t-1. This will be overwritten with the log-likelihood
* of the given emission probability up to time t.
* @param forwardLogProb Vector in which forward probabilities will be saved.
* Passing forwardLogProb as an empty vector indicates the start of the
* sequence (i.e. time t=0).
* @return Log-likelihood of the given sequence of emission up to time t.
*/
double EmissionLogLikelihood(const arma::vec& emissionLogProb,
double &logLikelihood,
arma::vec& forwardLogProb) const;
/**
* Compute the log of the scaling factor of the given data at time t.
* To calculate the log-likelihood for the whole sequence, accumulate the
* log scale factor (the return value of this function) over the entire
* sequence.
* This is meant for incremental or streaming computation of the
* log-likelihood of a sequence. For the first data point, provide an empty
* forwardLogProb vector.
*
* @param data observation at time t.
* @param forwardLogProb Vector in which forward probabilities will be saved.
* Passing forwardLogProb as an empty vector indicates the start of the
* sequence (i.e. time t=0).
* @return Log scale factor of the given sequence of data up at time t.
*/
double LogScaleFactor(const arma::vec &data,
arma::vec& forwardLogProb) const;
/**
* Compute the log-likelihood of the given data up to time t, storing the
* result in logLikelihood.
* This is meant for incremental or streaming computation of the
* log-likelihood of a sequence. For the first data point, provide an empty
* forwardLogProb vector.
*
* @param data observation at time t.
* @param logLikelihood Log-likelihood of the given sequence of data
* up to time t-1.
* @param forwardLogProb Vector in which forward probabilities will be saved.
* Passing forwardLogProb as an empty vector indicates the start of the
* sequence (i.e. time t=0).
* @return Log-likelihood of the given sequence of data up to time t.
*/
double LogLikelihood(const arma::vec &data,
double &logLikelihood,
arma::vec& forwardLogProb) const;
/**
* HMM filtering. Computes the k-step-ahead expected emission at each time
* conditioned only on prior observations. That is
@@ -366,6 +439,30 @@ class HMM
void save(Archive& ar, const uint32_t version) const;
protected:
/**
* Given emission probabilities, computes forward probabilities at time t=0.
*
* @param emissionLogProb Emission probability at time t=0.
* @param logScales Vector in which the log of scaling factors will be saved.
* @return Forward probabilities
*/
arma::vec ForwardAtT0(
const arma::vec& emissionLogProb,
double& logScales) const;
/**
* Given emission probabilities, computes forward probabilities for time t>0.
*
* @param emissionLogProb Emission probability at time t>0.
* @param logScales Vector in which the log of scaling factors will be saved.
* @param prevForwardLogProb Previous forward probabilities.
* @return Forward probabilities
*/
arma::vec ForwardAtTn(
const arma::vec& emissionLogProb,
double& logScales,
const arma::vec& prevForwardLogProb) const;
// Helper functions.
/**
* The Forward algorithm (part of the Forward-Backward algorithm). Computes
@@ -374,7 +471,7 @@ class HMM
* states and columns equal to the number of observations.
*
* @param dataSeq Data sequence to compute probabilities for.
* @param logScales Vector in which scaling factors will be saved.
* @param logScales Vector in which the log of scaling factors will be saved.
* @param forwardLogProb Matrix in which forward probabilities will be saved.
*/
void Forward(const arma::mat& dataSeq,
@@ -389,7 +486,7 @@ class HMM
* columns equal to the number of observations.
*
* @param dataSeq Data sequence to compute probabilities for.
* @param logScales Vector of scaling factors.
* @param logScales Vector of log of scaling factors.
* @param backwardLogProb Matrix in which backward probabilities will be saved.
*/
void Backward(const arma::mat& dataSeq,
+149 -26
View File
@@ -166,7 +166,7 @@ double HMM<Distribution>::Train(const std::vector<arma::mat>& dataSeq)
{
// Estimate of T_ij (probability of transition from state j to state
// i). We postpone multiplication of the old T_ij until later.
for (size_t i = 0; i < logTransition.n_rows; ++i)
for (size_t i = 0; i < logTransition.n_rows; i++)
{
newLogTransition(i, j) = math::LogAdd(newLogTransition(i, j),
forwardLog(j, t) + backwardLog(i, t + 1) +
@@ -204,7 +204,7 @@ double HMM<Distribution>::Train(const std::vector<arma::mat>& dataSeq)
logTransition += newLogTransition;
// Now we normalize the transition matrix.
for (size_t i = 0; i < logTransition.n_cols; ++i)
for (size_t i = 0; i < logTransition.n_cols; i++)
{
const double sum = math::AccuLog(logTransition.col(i));
if (std::isfinite(sum))
@@ -309,7 +309,7 @@ void HMM<Distribution>::Train(const std::vector<arma::mat>& dataSeq,
if (emissionList[state].size() > 0)
{
arma::mat emissions(dimensionality, emissionList[state].size());
for (size_t i = 0; i < emissions.n_cols; ++i)
for (size_t i = 0; i < emissions.n_cols; i++)
{
emissions.col(i) = dataSeq[emissionList[state][i].first].col(
emissionList[state][i].second);
@@ -486,7 +486,7 @@ double HMM<Distribution>::Predict(const arma::mat& dataSeq,
// Assemble the state probability for this element.
// Given that we are in state j, we use state with the highest probability
// of being the previous state.
for (size_t j = 0; j < logTransition.n_rows; ++j)
for (size_t j = 0; j < logTransition.n_rows; j++)
{
arma::vec prob = logStateProb.col(t - 1) + logTransition.row(j).t();
logStateProb(j, t) = prob.max(index) +
@@ -522,6 +522,80 @@ double HMM<Distribution>::LogLikelihood(const arma::mat& dataSeq) const
return accu(logScales);
}
/**
* Compute the log of the scaling factor of the given emission probability
* at time t. To calculate the log-likelihood for the whole sequence,
* accumulate log scale over the entire sequence
*/
template<typename Distribution>
double HMM<Distribution>::EmissionLogScaleFactor(
const arma::vec& emissionLogProb,
arma::vec& forwardLogProb) const
{
double curLogScale;
if (forwardLogProb.empty())
{
// We are at the start of the sequence (i.e. time t=0).
forwardLogProb = ForwardAtT0(emissionLogProb, curLogScale);
}
else
{
forwardLogProb = ForwardAtTn(emissionLogProb, curLogScale,
forwardLogProb);
}
return curLogScale;
}
/**
* Compute the log-likelihood of the given emission probability up to time t
*/
template<typename Distribution>
double HMM<Distribution>::EmissionLogLikelihood(
const arma::vec& emissionLogProb,
double& logLikelihood,
arma::vec& forwardLogProb) const
{
bool isStartOfSeq = forwardLogProb.empty();
double curLogScale = EmissionLogScaleFactor(emissionLogProb,
forwardLogProb);
logLikelihood = isStartOfSeq ? curLogScale : curLogScale + logLikelihood;
return logLikelihood;
}
/**
* Compute the log of the scaling factor of the given data at time t.
* To calculate the log-likelihood for the whole sequence, accumulate log
* scale over the entire sequence
*/
template<typename Distribution>
double HMM<Distribution>::LogScaleFactor(const arma::vec &data,
arma::vec& forwardLogProb) const
{
arma::vec emissionLogProb(logTransition.n_rows);
for (size_t state = 0; state < logTransition.n_rows; state++)
{
emissionLogProb(state) = emission[state].LogProbability(data);
}
return EmissionLogScaleFactor(emissionLogProb, forwardLogProb);
}
/**
* Compute the log-likelihood of the given data up to time t
*/
template<typename Distribution>
double HMM<Distribution>::LogLikelihood(const arma::vec& data,
double& logLikelihood,
arma::vec& forwardLogProb) const
{
bool isStartOfSeq = forwardLogProb.empty();
double curLogScale = LogScaleFactor(data, forwardLogProb);
logLikelihood = isStartOfSeq ? curLogScale : curLogScale + logLikelihood;
return logLikelihood;
}
/**
* HMM filtering.
*/
@@ -544,7 +618,7 @@ void HMM<Distribution>::Filter(const arma::mat& dataSeq,
// Compute expected emissions.
// Will not work for distributions without a Mean() function.
filterSeq.zeros(dimensionality, dataSeq.n_cols);
for (size_t i = 0; i < emission.size(); ++i)
for (size_t i = 0; i < emission.size(); i++)
filterSeq += emission[i].Mean() * forwardProb.row(i);
}
@@ -566,10 +640,67 @@ void HMM<Distribution>::Smooth(const arma::mat& dataSeq,
// Compute expected emissions.
// Will not work for distributions without a Mean() function.
smoothSeq.zeros(dimensionality, dataSeq.n_cols);
for (size_t i = 0; i < emission.size(); ++i)
for (size_t i = 0; i < emission.size(); i++)
smoothSeq += emission[i].Mean() * exp(stateLogProb.row(i));
}
/**
* The Forward procedure (part of the Forward-Backward algorithm).
*/
template<typename Distribution>
arma::vec HMM<Distribution>::ForwardAtT0(const arma::vec& emissionLogProb,
double& logScales) const
{
// Our goal is to calculate the forward probabilities:
// P(X_k | o_{1:k}) for all possible states X_k, for each time point k.
ConvertToLogSpace();
arma::vec forwardLogProb(logTransition.n_rows);
forwardLogProb.fill(-std::numeric_limits<double>::infinity());
// The first entry in the forward algorithm uses the initial state
// probabilities. Note that MATLAB assumes that the starting state (at
// t = -1) is state 0; this is not our assumption here. To force that
// behavior, you could append a single starting state to every single data
// sequence and that should produce results in line with MATLAB.
forwardLogProb = logInitial + emissionLogProb;
// Normalize probability.
logScales = math::AccuLog(forwardLogProb);
if (std::isfinite(logScales))
forwardLogProb -= logScales;
return forwardLogProb;
}
/**
* The Forward procedure (part of the Forward-Backward algorithm).
*/
template<typename Distribution>
arma::vec HMM<Distribution>::ForwardAtTn(const arma::vec& emissionLogProb,
double& logScales,
const arma::vec& prevForwardLogProb) const
{
// Our goal is to calculate the forward probabilities:
// P(X_k | o_{1:k}) for all possible states X_k, for each time point k.
arma::vec forwardLogProb(logTransition.n_rows);
forwardLogProb.fill(-std::numeric_limits<double>::infinity());
// Now compute the probabilities for each successive observation.
for (size_t state = 0; state < logTransition.n_rows; state++) {
// The forward probability of state j at time t is the sum over all
// states of the probability of the previous state transitioning to
// the current state and emitting the given observation.
arma::vec tmp = prevForwardLogProb + logTransition.row(state).t();
forwardLogProb(state) = math::AccuLog(tmp) + emissionLogProb(state);
}
// Normalize probability.
logScales = math::AccuLog(forwardLogProb);
if (std::isfinite(logScales))
forwardLogProb -= logScales;
return forwardLogProb;
}
/**
* The Forward procedure (part of the Forward-Backward algorithm).
*/
@@ -585,41 +716,32 @@ void HMM<Distribution>::Forward(const arma::mat& dataSeq,
logScales.resize(dataSeq.n_cols);
logScales.fill(-std::numeric_limits<double>::infinity());
ConvertToLogSpace();
// The first entry in the forward algorithm uses the initial state
// probabilities. Note that MATLAB assumes that the starting state (at
// t = -1) is state 0; this is not our assumption here. To force that
// behavior, you could append a single starting state to every single data
// sequence and that should produce results in line with MATLAB.
arma::vec emissionLogProb(logTransition.n_rows);
for (size_t state = 0; state < logTransition.n_rows; state++)
{
forwardLogProb(state, 0) = logInitial(state) +
emissionLogProb(state) =
emission[state].LogProbability(dataSeq.unsafe_col(0));
}
// Then normalize the column.
logScales[0] = math::AccuLog(forwardLogProb.col(0));
if (std::isfinite(logScales[0]))
forwardLogProb.col(0) -= logScales[0];
forwardLogProb.col(0) = ForwardAtT0(emissionLogProb, logScales(0));
// Now compute the probabilities for each successive observation.
for (size_t t = 1; t < dataSeq.n_cols; t++)
{
for (size_t j = 0; j < logTransition.n_rows; ++j)
for (size_t state = 0; state < logTransition.n_rows; state++)
{
// The forward probability of state j at time t is the sum over all states
// of the probability of the previous state transitioning to the current
// state and emitting the given observation.
arma::vec tmp = forwardLogProb.col(t - 1) + logTransition.row(j).t();
forwardLogProb(j, t) = math::AccuLog(tmp) +
emission[j].LogProbability(dataSeq.unsafe_col(t));
emissionLogProb(state) =
emission[state].LogProbability(dataSeq.unsafe_col(t));
}
// Normalize probability.
logScales[t] = math::AccuLog(forwardLogProb.col(t));
if (std::isfinite(logScales[t]))
forwardLogProb.col(t) -= logScales[t];
forwardLogProb.col(t) =
ForwardAtTn(emissionLogProb, logScales(t), forwardLogProb.col(t-1));
}
}
@@ -639,7 +761,7 @@ void HMM<Distribution>::Backward(const arma::mat& dataSeq,
// Now step backwards through all other observations.
for (size_t t = dataSeq.n_cols - 2; t + 1 > 0; t--)
{
for (size_t j = 0; j < logTransition.n_rows; ++j)
for (size_t j = 0; j < logTransition.n_rows; j++)
{
// The backward probability of state j at time t is the sum over all state
// of the probability of the next state having been a transition from the
@@ -660,7 +782,8 @@ void HMM<Distribution>::Backward(const arma::mat& dataSeq,
}
/**
* Make sure the variables in log space are in sync with the linear counter parts
* Make sure the variables in log space are in sync with the linear
* counterparts.
*/
template<typename Distribution>
void HMM<Distribution>::ConvertToLogSpace() const
+254 -3
View File
@@ -800,6 +800,7 @@ TEST_CASE("GaussianHMMPredictTest", "[HMMTest]")
emission.Covariance(cov.at(i));
}
// 100 2D observations.
arma::mat obs = {
{
-0.0424, -0.0395, -0.0336, -0.0294, -0.0299, -0.032, -0.0289, -0.0148,
@@ -832,9 +833,261 @@ TEST_CASE("GaussianHMMPredictTest", "[HMMTest]")
0.0521, 0.0313, 0.0188, 0.0113, 0.0068, 0.0042, 0.0026, 0.0018, 0.0014
}
};
//100 pre-calculated emission probabilities each for 10 states
std::vector<arma::vec> emissionProb = {
{ -2.7301e+03, 1.7874e+00, -1.9428e+00, -3.6365e+00, -4.0397e-01,
-1.5115e-01, -1.0328e+00, -1.1071e+00, 5.2876e-01, -1.0643e-01 },
{ -2.3684e+03, 1.8059e+00, -2.2058e+00, -4.0514e+00, -5.0935e-01,
-2.1126e-01, -1.1962e+00, -1.2567e+00, 4.1247e-01, -3.0199e-01 },
{ -1.7117e+03, 1.7981e+00, -2.5275e+00, -4.5634e+00, -6.4839e-01,
-2.9579e-01, -1.4000e+00, -1.4461e+00, 2.3795e-01, -5.5622e-01 },
{ -1.3089e+03, 1.7393e+00, -2.8685e+00, -5.0996e+00, -8.0288e-01,
-3.9863e-01, -1.6229e+00, -1.6478e+00, 2.3300e-02, -8.6617e-01 },
{ -1.3541e+03, 1.6414e+00, -3.1971e+00, -5.6043e+00, -9.5605e-01,
-5.1013e-01, -1.8460e+00, -1.8395e+00, -2.0603e-01, -1.2176e+00 },
{ -1.5521e+03, 1.5367e+00, -3.4806e+00, -6.0349e+00, -1.0924e+00,
-6.1426e-01, -2.0436e+00, -2.0051e+00, -4.2045e-01, -1.5500e+00 },
{ -1.2647e+03, 1.4680e+00, -3.6577e+00, -6.3144e+00, -1.1823e+00,
-6.8009e-01, -2.1646e+00, -2.1147e+00, -5.6512e-01, -1.7360e+00 },
{ -3.2650e+02, 1.4646e+00, -3.6693e+00, -6.3649e+00, -1.1957e+00,
-6.7711e-01, -2.1592e+00, -2.1377e+00, -5.8400e-01, -1.6543e+00 },
{ -1.3035e+02, 1.5123e+00, -3.5018e+00, -6.1593e+00, -1.1254e+00,
-6.0037e-01, -2.0181e+00, -2.0646e+00, -4.6413e-01, -1.3011e+00 },
{ -2.6279e+03, 1.5809e+00, -3.1559e+00, -5.6861e+00, -9.7699e-01,
-4.5903e-01, -1.7490e+00, -1.8956e+00, -2.2135e-01, -7.2772e-01 },
{ -9.6164e+03, 1.6193e+00, -2.6708e+00, -4.9944e+00, -7.8159e-01,
-2.8441e-01, -1.3909e+00, -1.6574e+00, 7.8411e-02, -5.0595e-02 },
{ -2.0944e+04, 1.5980e+00, -2.1094e+00, -4.1681e+00, -5.8143e-01,
-1.2105e-01, -1.0055e+00, -1.3879e+00, 3.4591e-01, 5.6464e-01 },
{ -3.3843e+04, 1.5241e+00, -1.5331e+00, -3.2977e+00, -4.1342e-01,
-8.0772e-03, -6.5026e-01, -1.1226e+00, 5.0244e-01, 9.8522e-01 },
{ -4.6678e+04, 1.3796e+00, -9.8223e-01, -2.4507e+00, -3.0368e-01,
3.0609e-02, -3.5787e-01, -8.8913e-01, 4.9234e-01, 1.1530e+00 },
{ -6.0839e+04, 1.1013e+00, -4.8302e-01, -1.6712e+00, -2.7541e-01,
-2.2814e-02, -1.4691e-01, -7.1095e-01, 2.6698e-01, 1.0338e+00 },
{ -7.8940e+04, 6.2341e-01, -6.4826e-02, -1.0034e+00, -3.5198e-01,
-1.8517e-01, -3.7803e-02, -6.1240e-01, -2.1704e-01, 5.8353e-01 },
{ -1.0182e+05, -8.9362e-02, 2.5888e-01, -4.6429e-01, -5.5297e-01,
-4.7752e-01, -4.9871e-02, -6.0739e-01, -1.0089e+00, -2.6587e-01 },
{ -1.2437e+05, -9.8625e-01, 4.7236e-01, -8.1256e-02, -8.8097e-01,
-9.0979e-01, -2.0039e-01, -6.9837e-01, -2.1229e+00, -1.5424e+00 },
{ -1.3878e+05, -1.9546e+00, 5.6976e-01, 1.2361e-01, -1.3043e+00,
-1.4534e+00, -4.7690e-01, -8.6807e-01, -3.4831e+00, -3.1393e+00 },
{ -1.3979e+05, -2.8896e+00, 5.6962e-01, 1.6577e-01, -1.7631e+00,
-2.0456e+00, -8.3102e-01, -1.0792e+00, -4.9380e+00, -4.8430e+00 },
{ -1.2717e+05, -3.7474e+00, 5.0493e-01, 9.2444e-02, -2.1969e+00,
-2.6201e+00, -1.2028e+00, -1.2907e+00, -6.3319e+00, -6.4416e+00 },
{ -1.0548e+05, -4.5565e+00, 4.0397e-01, -4.5775e-02, -2.5711e+00,
-3.1354e+00, -1.5493e+00, -1.4771e+00, -7.5697e+00, -7.8170e+00 },
{ -8.0621e+04, -5.3691e+00, 2.8365e-01, -2.1252e-01, -2.8783e+00,
-3.5784e+00, -1.8523e+00, -1.6312e+00, -8.6245e+00, -8.9480e+00 },
{ -5.6310e+04, -6.2411e+00, 1.5008e-01, -3.9022e-01, -3.1294e+00,
-3.9597e+00, -2.1142e+00, -1.7569e+00, -9.5239e+00, -9.8785e+00 },
{ -3.4173e+04, -7.2306e+00, 3.0396e-03, -5.7242e-01, -3.3347e+00,
-4.2928e+00, -2.3417e+00, -1.8583e+00, -1.0301e+01, -1.0652e+01 },
{ -1.5877e+04, -8.3900e+00, -1.5871e-01, -7.5362e-01, -3.4959e+00,
-4.5816e+00, -2.5356e+00, -1.9353e+00, -1.0963e+01, -1.1284e+01 },
{ -3.3829e+03, -9.7572e+00, -3.3006e-01, -9.1554e-01, -3.5912e+00,
-4.8035e+00, -2.6770e+00, -1.9722e+00, -1.1452e+01, -1.1714e+01 },
{ -5.6088e+02, -1.1394e+01, -5.0305e-01, -1.0261e+00, -3.5777e+00,
-4.9138e+00, -2.7301e+00, -1.9403e+00, -1.1653e+01, -1.1829e+01 },
{ -1.4303e+04, -1.3346e+01, -6.7336e-01, -1.0564e+00, -3.4266e+00,
-4.8757e+00, -2.6690e+00, -1.8219e+00, -1.1470e+01, -1.1561e+01 },
{ -4.9066e+04, -1.5534e+01, -8.4176e-01, -1.0079e+00, -3.1636e+00,
-4.7028e+00, -2.5116e+00, -1.6369e+00, -1.0937e+01, -1.0995e+01 },
{ -9.9717e+04, -1.7702e+01, -1.0039e+00, -9.1443e-01, -2.8597e+00,
-4.4595e+00, -2.3138e+00, -1.4339e+00, -1.0224e+01, -1.0331e+01 },
{ -1.5886e+05, -1.9676e+01, -1.1535e+00, -7.9762e-01, -2.5479e+00,
-4.1805e+00, -2.1039e+00, -1.2332e+00, -9.4233e+00, -9.6530e+00 },
{ -2.2947e+05, -2.1635e+01, -1.3117e+00, -6.6325e-01, -2.2133e+00,
-3.8587e+00, -1.8780e+00, -1.0253e+00, -8.5051e+00, -8.9416e+00 },
{ -3.1968e+05, -2.3792e+01, -1.5095e+00, -5.1672e-01, -1.8381e+00,
-3.4770e+00, -1.6312e+00, -8.0190e-01, -7.4108e+00, -8.1836e+00 },
{ -4.3323e+05, -2.6183e+01, -1.7728e+00, -3.8390e-01, -1.4394e+00,
-3.0487e+00, -1.3857e+00, -5.7953e-01, -6.1647e+00, -7.4521e+00 },
{ -5.6473e+05, -2.8589e+01, -2.1061e+00, -3.0168e-01, -1.0547e+00,
-2.6054e+00, -1.1773e+00, -3.8708e-01, -4.8475e+00, -6.8476e+00 },
{ -6.9974e+05, -3.0612e+01, -2.4921e+00, -3.0913e-01, -7.2535e-01,
-2.1849e+00, -1.0419e+00, -2.5359e-01, -3.5677e+00, -6.4479e+00 },
{ -8.0655e+05, -3.1539e+01, -2.8524e+00, -4.2185e-01, -4.8692e-01,
-1.8260e+00, -9.9484e-01, -1.9514e-01, -2.4629e+00, -6.2373e+00 },
{ -8.5216e+05, -3.0655e+01, -3.0833e+00, -6.1881e-01, -3.3169e-01,
-1.5249e+00, -1.0091e+00, -1.9595e-01, -1.5717e+00, -6.0513e+00 },
{ -8.2392e+05, -2.7811e+01, -3.1362e+00, -8.7526e-01, -2.3459e-01,
-1.2631e+00, -1.0480e+00, -2.3278e-01, -8.7344e-01, -5.7341e+00 },
{ -7.3612e+05, -2.3582e+01, -3.0425e+00, -1.1744e+00, -1.7841e-01,
-1.0351e+00, -1.0893e+00, -2.9210e-01, -3.4780e-01, -5.2495e+00 },
{ -6.1397e+05, -1.8706e+01, -2.8744e+00, -1.5195e+00, -1.5516e-01,
-8.3816e-01, -1.1304e+00, -3.7330e-01, 3.7744e-02, -4.6424e+00 },
{ -4.8041e+05, -1.3799e+01, -2.7054e+00, -1.9262e+00, -1.6637e-01,
-6.7558e-01, -1.1810e+00, -4.8405e-01, 3.0197e-01, -3.9898e+00 },
{ -3.4790e+05, -9.2300e+00, -2.5518e+00, -2.3683e+00, -2.0524e-01,
-5.4582e-01, -1.2297e+00, -6.1666e-01, 4.5657e-01, -3.3063e+00 },
{ -2.2370e+05, -5.1887e+00, -2.3941e+00, -2.7911e+00, -2.5500e-01,
-4.3560e-01, -1.2487e+00, -7.5224e-01, 5.3161e-01, -2.5570e+00 },
{ -1.1273e+05, -1.7195e+00, -2.2258e+00, -3.1794e+00, -3.0915e-01,
-3.2974e-01, -1.2221e+00, -8.8867e-01, 5.5755e-01, -1.7017e+00 },
{ -2.8363e+04, 9.0588e-01, -2.0601e+00, -3.5233e+00, -3.7171e-01,
-2.2162e-01, -1.1370e+00, -1.0334e+00, 5.4434e-01, -7.3209e-01 },
{ -1.2122e+03, 1.9784e+00, -1.9455e+00, -3.7971e+00, -4.5862e-01,
-1.2081e-01, -9.8979e-01, -1.2000e+00, 4.7839e-01, 2.5783e-01 },
{ -7.1694e+04, 5.6327e-01, -2.0051e+00, -4.0345e+00, -6.1306e-01,
-6.6602e-02, -8.2833e-01, -1.4287e+00, 3.0684e-01, 1.0022e+00 },
{ -2.6198e+05, -3.8345e+00, -2.3396e+00, -4.2798e+00, -8.6900e-01,
-9.8822e-02, -7.1489e-01, -1.7486e+00, -1.6219e-02, 1.2358e+00 },
{ -5.5328e+05, -1.0687e+01, -2.9124e+00, -4.5058e+00, -1.2121e+00,
-2.2259e-01, -6.7273e-01, -2.1347e+00, -4.7585e-01, 8.8080e-01 },
{ -8.9436e+05, -1.8602e+01, -3.5518e+00, -4.6037e+00, -1.5911e+00,
-4.1140e-01, -6.7958e-01, -2.5173e+00, -1.0137e+00, 3.7886e-02 },
{ -1.2162e+06, -2.5781e+01, -4.0541e+00, -4.4848e+00, -1.9485e+00,
-6.3137e-01, -7.0903e-01, -2.8240e+00, -1.5699e+00, -1.1063e+00 },
{ -1.4436e+06, -3.0414e+01, -4.2395e+00, -4.1197e+00, -2.2265e+00,
-8.4654e-01, -7.3852e-01, -2.9921e+00, -2.0869e+00, -2.2970e+00 },
{ -1.5227e+06, -3.1337e+01, -3.9989e+00, -3.5197e+00, -2.3836e+00,
-1.0315e+00, -7.4887e-01, -2.9823e+00, -2.5313e+00, -3.3017e+00 },
{ -1.4386e+06, -2.8472e+01, -3.3472e+00, -2.7563e+00, -2.4087e+00,
-1.1801e+00, -7.3524e-01, -2.7971e+00, -2.9034e+00, -3.9803e+00 },
{ -1.2257e+06, -2.2958e+01, -2.4364e+00, -1.9521e+00, -2.3275e+00,
-1.3043e+00, -7.0875e-01, -2.4858e+00, -3.2295e+00, -4.3252e+00 },
{ -9.4813e+05, -1.6527e+01, -1.4675e+00, -1.2121e+00, -2.1965e+00,
-1.4367e+00, -6.9389e-01, -2.1228e+00, -3.5740e+00, -4.4844e+00 },
{ -6.6589e+05, -1.0680e+01, -6.1313e-01, -6.1440e-01, -2.0638e+00,
-1.5984e+00, -7.0917e-01, -1.7727e+00, -3.9726e+00, -4.5979e+00 },
{ -4.1809e+05, -6.2975e+00, 3.1651e-02, -1.8731e-01, -1.9586e+00,
-1.7982e+00, -7.6241e-01, -1.4730e+00, -4.4365e+00, -4.7645e+00 },
{ -2.2534e+05, -3.7546e+00, 4.3188e-01, 7.1872e-02, -1.8959e+00,
-2.0366e+00, -8.5455e-01, -1.2417e+00, -4.9637e+00, -5.0424e+00 },
{ -9.4330e+04, -3.0403e+00, 5.9517e-01, 1.8422e-01, -1.8702e+00,
-2.2952e+00, -9.7314e-01, -1.0776e+00, -5.5109e+00, -5.4155e+00 },
{ -2.1454e+04, -3.9202e+00, 5.5647e-01, 1.8381e-01, -1.8704e+00,
-2.5578e+00, -1.1056e+00, -9.6899e-01, -6.0419e+00, -5.8579e+00 },
{ -31.4830, -6.0953, 0.3567, 0.1044, -1.8840, -2.8086, -1.2397,
-0.9026, -6.5224, -6.3374 },
{ -2.2442e+04, -9.2735e+00, 3.4960e-02, -2.1605e-02, -1.8931e+00,
-3.0282e+00, -1.3611e+00, -8.6066e-01, -6.9076e+00, -6.8075e+00 },
{ -8.1676e+04, -1.3138e+01, -3.6831e-01, -1.6104e-01, -1.8763e+00,
-3.1905e+00, -1.4522e+00, -8.2511e-01, -7.1362e+00, -7.2081e+00 },
{ -1.6865e+05, -1.7287e+01, -8.0643e-01, -2.8264e-01, -1.8178e+00,
-3.2726e+00, -1.4987e+00, -7.8144e-01, -7.1585e+00, -7.4877e+00 },
{ -2.7001e+05, -2.1213e+01, -1.2247e+00, -3.6116e-01, -1.7095e+00,
-3.2596e+00, -1.4928e+00, -7.2002e-01, -6.9485e+00, -7.6058e+00 },
{ -3.7506e+05, -2.4628e+01, -1.5962e+00, -3.9394e-01, -1.5583e+00,
-3.1610e+00, -1.4428e+00, -6.4101e-01, -6.5350e+00, -7.5763e+00 },
{ -4.7871e+05, -2.7455e+01, -1.9194e+00, -3.9090e-01, -1.3720e+00,
-2.9900e+00, -1.3606e+00, -5.4763e-01, -5.9492e+00, -7.4279e+00 },
{ -5.7329e+05, -2.9501e+01, -2.1830e+00, -3.6323e-01, -1.1594e+00,
-2.7564e+00, -1.2564e+00, -4.4501e-01, -5.2194e+00, -7.1738e+00 },
{ -6.4968e+05, -3.0560e+01, -2.3747e+00, -3.2775e-01, -9.3375e-01,
-2.4742e+00, -1.1428e+00, -3.4141e-01, -4.3880e+00, -6.8281e+00 },
{ -6.9933e+05, -3.0501e+01, -2.4875e+00, -3.0631e-01, -7.1262e-01,
-2.1653e+00, -1.0343e+00, -2.4789e-01, -3.5174e+00, -6.4120e+00 },
{ -7.1802e+05, -2.9350e+01, -2.5271e+00, -3.2061e-01, -5.1194e-01,
-1.8521e+00, -9.4328e-01, -1.7450e-01, -2.6686e+00, -5.9486e+00 },
{ -7.0553e+05, -2.7236e+01, -2.5060e+00, -3.8730e-01, -3.4217e-01,
-1.5515e+00, -8.7707e-01, -1.2819e-01, -1.8857e+00, -5.4542e+00 },
{ -6.6569e+05, -2.4393e+01, -2.4435e+00, -5.1663e-01, -2.1031e-01,
-1.2775e+00, -8.3941e-01, -1.1339e-01, -1.2023e+00, -4.9470e+00 },
{ -6.1301e+05, -2.1269e+01, -2.3992e+00, -7.3370e-01, -1.2064e-01,
-1.0383e+00, -8.4300e-01, -1.3855e-01, -6.1878e-01, -4.4864e+00 },
{ -5.6195e+05, -1.8233e+01, -2.4507e+00, -1.0921e+00, -8.5743e-02,
-8.4467e-01, -9.1749e-01, -2.2378e-01, -1.3441e-01, -4.1677e+00 },
{ -5.0308e+05, -1.5078e+01, -2.5824e+00, -1.6122e+00, -1.1850e-01,
-7.0720e-01, -1.0690e+00, -3.7916e-01, 2.0948e-01, -3.9737e+00 },
{ -4.2417e+05, -1.1613e+01, -2.7333e+00, -2.2592e+00, -2.1392e-01,
-6.2529e-01, -1.2691e+00, -5.9232e-01, 3.8196e-01, -3.8021e+00 },
{ -3.4311e+05, -8.4490e+00, -2.9262e+00, -2.9840e+00, -3.6201e-01,
-6.0612e-01, -1.5036e+00, -8.4657e-01, 3.8172e-01, -3.6994e+00 },
{ -2.6553e+05, -5.7959e+00, -3.0657e+00, -3.6135e+00, -5.0450e-01,
-6.1056e-01, -1.6893e+00, -1.0726e+00, 2.9263e-01, -3.5310e+00 },
{ -1.6581e+05, -2.8806e+00, -2.9242e+00, -3.9480e+00, -5.5108e-01,
-5.3603e-01, -1.6743e+00, -1.1832e+00, 2.8121e-01, -2.8215e+00 },
{ -6.3112e+04, -4.4355e-02, -2.4673e+00, -3.8848e+00, -4.8010e-01,
-3.5415e-01, -1.4075e+00, -1.1547e+00, 4.0803e-01, -1.5227e+00 },
{ -5.4196e+03, 1.6750e+00, -1.9272e+00, -3.5655e+00, -3.8433e-01,
-1.5578e-01, -1.0312e+00, -1.0745e+00, 5.4628e-01, -1.8838e-01 },
{ -7.9742e+03, 1.9542e+00, -1.5297e+00, -3.2224e+00, -3.5023e-01,
-2.9557e-02, -7.1541e-01, -1.0340e+00, 5.7335e-01, 7.0234e-01 },
{ -4.6838e+04, 1.3383e+00, -1.2840e+00, -2.9202e+00, -3.6943e-01,
2.1035e-02, -4.9879e-01, -1.0295e+00, 5.0388e-01, 1.1296e+00 },
{ -9.2965e+04, 5.0293e-01, -1.1033e+00, -2.6251e+00, -4.0461e-01,
2.4992e-02, -3.5062e-01, -1.0246e+00, 3.8909e-01, 1.2595e+00 },
{ -1.3250e+05, -2.3738e-01, -9.9398e-01, -2.4136e+00, -4.4740e-01,
6.0968e-03, -2.6559e-01, -1.0308e+00, 2.6684e-01, 1.2440e+00 },
{ -1.7149e+05, -9.7999e-01, -9.1698e-01, -2.2384e+00, -4.9912e-01,
-2.5475e-02, -2.0762e-01, -1.0468e+00, 1.3078e-01, 1.1599e+00 },
{ -2.2091e+05, -1.9350e+00, -8.5497e-01, -2.0582e+00, -5.7508e-01,
-7.7816e-02, -1.6138e-01, -1.0794e+00, -5.6045e-02, 9.8875e-01 },
{ -2.8140e+05, -3.1219e+00, -8.2568e-01, -1.8962e+00, -6.7862e-01,
-1.5242e-01, -1.3554e-01, -1.1353e+00, -2.9320e-01, 7.2082e-01 },
{ -3.4167e+05, -4.3171e+00, -8.2824e-01, -1.7733e+00, -7.8907e-01,
-2.3483e-01, -1.3167e-01, -1.2015e+00, -5.3627e-01, 4.0854e-01 },
{ -3.7868e+05, -5.0537e+00, -8.3691e-01, -1.7046e+00, -8.6035e-01,
-2.9036e-01, -1.3690e-01, -1.2447e+00, -6.9304e-01, 1.9260e-01 },
{ -3.7429e+05, -4.9406e+00, -7.8456e-01, -1.6323e+00, -8.6203e-01,
-3.0644e-01, -1.3159e-01, -1.2267e+00, -7.3358e-01, 1.3468e-01 },
{ -3.3293e+05, -4.0758e+00, -6.6873e-01, -1.5416e+00, -8.0032e-01,
-2.8877e-01, -1.1346e-01, -1.1498e+00, -6.7458e-01, 2.1365e-01 },
{ -2.7541e+05, -2.9085e+00, -5.3210e-01, -1.4470e+00, -7.0706e-01,
-2.5517e-01, -9.1435e-02, -1.0445e+00, -5.6402e-01, 3.5113e-01 },
{ -2.2010e+05, -1.8209e+00, -4.1116e-01, -1.3627e+00, -6.1220e-01,
-2.2144e-01, -7.3319e-02, -9.4005e-01, -4.4660e-01, 4.7992e-01 },
{ -1.7809e+05, -1.0242e+00, -3.2646e-01, -1.3011e+00, -5.3612e-01,
-1.9567e-01, -6.2291e-02, -8.5731e-01, -3.5032e-01, 5.7022e-01 },
{ -1.5426e+05, -5.8691e-01, -2.8121e-01, -1.2660e+00, -4.9111e-01,
-1.8141e-01, -5.7387e-02, -8.0842e-01, -2.9317e-01, 6.1601e-01 },
};
const double loglikelihoodRef = -2734.43;
// Test log-likelihood calculation for the whole data.
{
const double loglikelihood = hmm.LogLikelihood(obs);
REQUIRE(loglikelihood == Approx(loglikelihoodRef).epsilon(1e-3));
}
// Test loglikelihood calculation in an incremental way.
// It simulates the case where we have a stream of data.
{
double loglikelihood;
arma::vec forwardLogProb;
for (size_t t = 0; t<obs.n_cols; ++t)
{
loglikelihood = hmm.LogLikelihood(obs.col(t), loglikelihood,
forwardLogProb);
}
REQUIRE(loglikelihood == Approx(loglikelihoodRef).epsilon(1e-3));
}
// Test loglikelihood calculation in an incremental way.
// It simulates the case where we have a stream of data.
// In this case the accumulation of the log scales factor to calculate
// the log-likelihood value is done outside of the loop
{
double loglikelihood = 0;
arma::vec forwardLogProb;
for (size_t t = 0; t<obs.n_cols; ++t)
{
double logScale = hmm.LogScaleFactor(obs.col(t), forwardLogProb);
loglikelihood += logScale;
}
REQUIRE(loglikelihood == Approx(loglikelihoodRef).epsilon(1e-3));
}
// Test loglikelihood calculation in an incremental way.
// It simulates the case where we have emission probabilities pre-calculated.
{
double loglikelihood = 0;
arma::vec forwardLogProb;
for (size_t t = 0; t<emissionProb.size(); ++t)
{
loglikelihood = hmm.EmissionLogLikelihood(emissionProb.at(t),
loglikelihood, forwardLogProb);
}
REQUIRE(loglikelihood == Approx(loglikelihoodRef).epsilon(1e-1));
}
arma::Row<size_t> stateSeq;
auto likelihood = hmm.LogLikelihood(obs);
hmm.Predict(obs, stateSeq);
arma::Row<size_t> stateSeqRef = { 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
@@ -843,8 +1096,6 @@ TEST_CASE("GaussianHMMPredictTest", "[HMMTest]")
7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9,
9, 9, 9, 9, 9, 9, 9, 9, 9, 9 };
REQUIRE(likelihood == Approx(-2734.43).epsilon(1e-5));
for (size_t i = 0; i < stateSeqRef.n_cols; ++i)
{
REQUIRE(stateSeqRef.at(i) == stateSeq.at(i));