code: private double[][][] estimateXi(int[] o, double[][] alpha, double[][] beta) {
          if (o.length <= 1) {
              throw new IllegalArgumentException("Observation sequence is too short.");
          }
  
          int n = o.length - 1;
          double xi[][][] = new double[n][numStates][numStates];
  
          for (int t = 0; t < n; t++) {
              for (int i = 0; i < numStates; i++) {
                  for (int j = 0; j < numStates; j++) {
                      xi[t][i][j] = alpha[t][i] * a[i][j] * b[j][o[t + 1]] * beta[t + 1][j];
                  }
              }
          }
  
          return xi;
      }
docs: Here, the xi (and, thus, gamma) values are not divided by the probability
  of the sequence because this probability might be too small and induce an
  underflow. xi[t][i][j] still can be interpreted as P(q_t = i and q_(t+1)
  = j | O, HMM) because we assume that the scaling factors are such that
  their product is equal to the inverse of the probability of the sequence.