diff --git a/HISTORY.md b/HISTORY.md index 7529a3ee86..c8a3c94497 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -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 diff --git a/src/mlpack/methods/hmm/hmm.hpp b/src/mlpack/methods/hmm/hmm.hpp index eceb51b527..1821083a0e 100644 --- a/src/mlpack/methods/hmm/hmm.hpp +++ b/src/mlpack/methods/hmm/hmm.hpp @@ -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, diff --git a/src/mlpack/methods/hmm/hmm_impl.hpp b/src/mlpack/methods/hmm/hmm_impl.hpp index 8e4d8a2b2f..d3d1d72e61 100644 --- a/src/mlpack/methods/hmm/hmm_impl.hpp +++ b/src/mlpack/methods/hmm/hmm_impl.hpp @@ -166,7 +166,7 @@ double HMM::Train(const std::vector& 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::Train(const std::vector& 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::Train(const std::vector& 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::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::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 +double HMM::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 +double HMM::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 +double HMM::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 +double HMM::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::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::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 +arma::vec HMM::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::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 +arma::vec HMM::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::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::Forward(const arma::mat& dataSeq, logScales.resize(dataSeq.n_cols); logScales.fill(-std::numeric_limits::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::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::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 void HMM::ConvertToLogSpace() const diff --git a/src/mlpack/tests/hmm_test.cpp b/src/mlpack/tests/hmm_test.cpp index da1debd3a2..03a22f31e6 100644 --- a/src/mlpack/tests/hmm_test.cpp +++ b/src/mlpack/tests/hmm_test.cpp @@ -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 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 stateSeq; - auto likelihood = hmm.LogLikelihood(obs); hmm.Predict(obs, stateSeq); arma::Row 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));