public static HMM fit(int[][] observations, int[][] labels) {
        if (observations.length != labels.length) {
            throw new IllegalArgumentException("The number of observation sequences and that of label sequences are different.");
        }

        int N = 0; // the number of states
        int M = 0; // the number of symbols

        for (int i = 0; i < observations.length; i++) {
            if (observations[i].length != labels[i].length) {
                throw new IllegalArgumentException(String.format("The length of observation sequence %d and that of corresponding label sequence are different.", i));
            }

            N = Math.max(N, MathEx.max(labels[i]) + 1);
            M = Math.max(M, MathEx.max(observations[i]) + 1);
        }

        double[] pi = new double[N];
        double[][] a = new double[N][N];
        double[][] b = new double[N][M];

        for (int i = 0; i < observations.length; i++) {
            pi[labels[i][0]]++;
            b[labels[i][0]][observations[i][0]]++;
            for (int j = 1; j < observations[i].length; j++) {
                a[labels[i][j - 1]][labels[i][j]]++;
                b[labels[i][j]][observations[i][j]]++;
            }
        }

        MathEx.unitize1(pi);
        for (int i = 0; i < N; i++) {
            MathEx.unitize1(a[i]);
            MathEx.unitize1(b[i]);
        }

        return new HMM(pi, Matrix.of(a), Matrix.of(b));
    }
@Test
    public void testFit() {
        System.out.println("fit");
        MathEx.setSeed(19650218); // to get repeatable results.

        EmpiricalDistribution initial = new EmpiricalDistribution(pi);

        EmpiricalDistribution[] transition = new EmpiricalDistribution[a.length];
        for (int i = 0; i < transition.length; i++) {
            transition[i] = new EmpiricalDistribution(a[i]);
        }

        EmpiricalDistribution[] emission = new EmpiricalDistribution[b.length];
        for (int i = 0; i < emission.length; i++) {
            emission[i] = new EmpiricalDistribution(b[i]);
        }

        int[][] sequences = new int[5000][];
        int[][] labels = new int[5000][];
        for (int i = 0; i < sequences.length; i++) {
            sequences[i] = new int[30 * (MathEx.randomInt(5) + 1)];
            labels[i] = new int[sequences[i].length];
            int state = (int) initial.rand();
            sequences[i][0] = (int) emission[state].rand();
            labels[i][0] = state;
            for (int j = 1; j < sequences[i].length; j++) {
                state = (int) transition[state].rand();
                sequences[i][j] = (int) emission[state].rand();
                labels[i][j] = state;
            }
        }

        HMM model = HMM.fit(sequences, labels);
        System.out.println(model);

        double[] expPi2 = {0.5076, 0.4924};
        double[][] expA2 = {{0.8002, 0.1998}, {0.1987, 0.8013}};
        double[][] expB2 = {{0.5998, 0.4002}, {0.4003, 0.5997}};

        double[] pi2 = model.getInitialStateProbabilities();
        for (int i = 0; i < pi.length; i++) {
            assertEquals(expPi2[i], pi2[i], 1E-4);
        }

        Matrix a2 = model.getStateTransitionProbabilities();
        for (int i = 0; i < a.length; i++) {
            for (int j = 0; j < a[i].length; j++) {
                assertEquals(expA2[i][j], a2.get(i, j), 1E-4);
            }
        }

        Matrix b2 = model.getSymbolEmissionProbabilities();
        for (int i = 0; i < b.length; i++) {
            for (int j = 0; j < b[i].length; j++) {
                assertEquals(expB2[i][j], b2.get(i, j), 1E-4);
            }
        }
    }