/// <summary>
        ///   Predicts the next observation occurring after a given observation sequence.
        /// </summary>
        public double Predict(double[] observations)
        {
            double probability;
            double[][] obs = convert(observations);
            return Predict(obs, 1, out probability)[0][0];
        }

        /// <summary>
        ///   Predicts the next observation occurring after a given observation sequence.
        /// </summary>
        public double[] Predict(double[][] observations)
        {
            double probability;
            return Predict(observations, 1, out probability)[0];
        }

        /// <summary>
        ///   Predicts the next observation occurring after a given observation sequence.
        /// </summary>
        public double[] Predict(double[] observations, int next)
        {
            double probability;
            double[][] obs = convert(observations);
            return Predict(obs, next, out probability)[0];
        }

        /// <summary>
        ///   Predicts the next observation occurring after a given observation sequence.
        /// </summary>
        public double[][] Predict(double[][] observations, int next)
        {
            double probability;
            return Predict(observations, next, out probability);
        }

        /// <summary>
        ///   Predicts the next observation occurring after a given observation sequence.
        /// </summary>
        public double[] Predict(double[] observations, int next, out double probability)
        {
            double[][] obs = convert(observations);
            double[][] prediction = Predict(obs, next, out probability);
            return Accord.Math.Matrix.Combine(prediction);
        }

        /// <summary>
        ///   Predicts the next observation occurring after a given observation sequence.
        /// </summary>
        public double[][] Predict(double[][] observations, int next, out double probability)
        {
            int states = States;
            int T = next;

            double[,] A = Transitions;
            double logLikelihood;

            double[][] prediction = new double[next][];
            double[][] means = new double[states][];

            if (Dimension > 1)
            {
                for (int i = 0; i < states; i++)
                    means[i] = (B[i] as IMultivariateDistribution).Mean;
            }
            else
            {
                for (int i = 0; i < states; i++)
                    means[i] = new double[] { (B[i] as IUnivariateDistribution).Mean };
            }


            // Compute forward probabilities for the given observation sequence.
            double[,] fw0 = ForwardBackwardAlgorithm.Forward(this, observations, out logLikelihood);

            // Create a matrix to store the future probabilities for the prediction
            // sequence and copy the latest forward probabilities on its first row.
            double[,] fwd = new double[T + 1, States];


            // 1. Initialization
            for (int i = 0; i < States; i++)
                fwd[0, i] = fw0[observations.Length - 1, i];


            // 2. Induction
            for (int t = 0; t < T; t++)
            {
                double scale = 0;
                for (int i = 0; i < states; i++)
                {
                    double sum = 0.0;
                    for (int j = 0; j < states; j++)
                        sum += fwd[t, j] * A[j, i];
                    fwd[t + 1, i] = sum * B[i].ProbabilityFunction(means[i]);

                    scale += fwd[t + 1, i];
                }

                scale /= Math.Exp(logLikelihood);

                if (scale != 0) // Scaling
                {
                    for (int i = 0; i < states; i++)
                        fwd[t + 1, i] /= scale;
                }


                // Select most probable value
                double maxWeight = fwd[t + 1, 0];
                double sumWeight = maxWeight;
                prediction[t] = means[0];
                for (int i = 1; i < states; i++)
                {
                    if (fwd[t + 1, i] > maxWeight)
                    {
                        maxWeight = fwd[t + 1, i];
                        prediction[t] = means[i];
                    }
                }

                // Recompute log-likelihood
                logLikelihood = Math.Log(maxWeight);
            }

            // Returns the sequence probability as an out parameter
            probability = Math.Exp(logLikelihood);

            return prediction;
        }