for (int i = 0; i < numHiddenStates; i++) {
hiddenPriorProbability.set(i, random.nextDouble());
}
normalizeProbabilities();
DoubleMatrix alpha = new DenseDoubleMatrix(features.length, numHiddenStates);
DoubleMatrix beta = new DenseDoubleMatrix(features.length, numHiddenStates);
for (int iteration = 0; iteration < maxIterations; iteration++) {
// in every iteration we initialize a new model that is a copy of the old
DoubleMatrix transitionProbabilityMatrix = this.transitionProbabilityMatrix
.deepCopy();
DoubleMatrix emissionProbabilityMatrix = this.emissionProbabilityMatrix
.deepCopy();
DoubleVector hiddenPriorProbability = this.hiddenPriorProbability
.deepCopy();
// expectation step
alpha = forward(alpha, transitionProbabilityMatrix,
emissionProbabilityMatrix, hiddenPriorProbability, features);
beta = backward(beta, transitionProbabilityMatrix,
emissionProbabilityMatrix, hiddenPriorProbability, features);
// now do the real baum-welch algorithm / maximization step calculate the
// prior out of the alpha and beta factors in their first row
hiddenPriorProbability = alpha.getRowVector(0).multiply(
beta.getRowVector(0));
final double modelLikelihood = estimateLikelihood(alpha);
// compute real transition probabilities
for (int i = 0; i < numHiddenStates; i++) {
for (int j = 0; j < numHiddenStates; j++) {
double temp = 0d;
for (int t = 0; t < features.length - 1; t++) {
Iterator<DoubleVectorElement> iterateNonZero = features[t + 1]
.iterateNonZero();
while (iterateNonZero.hasNext()) {
temp += alpha.get(t, i)
* emissionProbabilityMatrix.get(j, iterateNonZero.next()
.getIndex()) * beta.get(t + 1, j);
}
}
transitionProbabilityMatrix.set(i, j,
transitionProbabilityMatrix.get(i, j) * temp / modelLikelihood);
}
}
// compute real emission probabilities
for (int i = 0; i < numHiddenStates; i++) {
for (int j = 0; j < numVisibleStates; j++) {
double temp = 0d;
for (int t = 0; t < features.length; t++) {
Iterator<DoubleVectorElement> iterateNonZero = features[t]
.iterateNonZero();
while (iterateNonZero.hasNext()) {
DoubleVectorElement next = iterateNonZero.next();
if (next.getIndex() == j) {
temp += alpha.get(t, i) * beta.get(t, i);
}
}
}
emissionProbabilityMatrix.set(i, j, temp / modelLikelihood);
}
}
// normalize again after baumwelch pass
normalize(hiddenPriorProbability, transitionProbabilityMatrix,