This function computes the statistics used to train an HMM joint segmentation-classification model
    using a sequence of sequential features and respective labels

    ARGUMENTS:
     - features:    a numpy matrix of feature vectors (numOfDimensions x n_wins)
     - labels:    a numpy array of class indices (n_wins x 1)
    RETURNS:
     - start_prob:    matrix of prior class probabilities (n_classes x 1)
     - transmat:    transition matrix (n_classes x n_classes)
     - means:    means matrix (numOfDimensions x 1)
     - cov:        deviation matrix (numOfDimensions x 1)

def trainHMM_computeStatistics(features, labels):
    '''
    This function computes the statistics used to train an HMM joint segmentation-classification model
    using a sequence of sequential features and respective labels

    ARGUMENTS:
     - features:    a numpy matrix of feature vectors (numOfDimensions x n_wins)
     - labels:    a numpy array of class indices (n_wins x 1)
    RETURNS:
     - start_prob:    matrix of prior class probabilities (n_classes x 1)
     - transmat:    transition matrix (n_classes x n_classes)
     - means:    means matrix (numOfDimensions x 1)
     - cov:        deviation matrix (numOfDimensions x 1)
    '''
    u_labels = numpy.unique(labels)
    n_comps = len(u_labels)

    n_feats = features.shape[0]

    if features.shape[1] < labels.shape[0]:
        print("trainHMM warning: number of short-term feature vectors "
              "must be greater or equal to the labels length!")
        labels = labels[0:features.shape[1]]

    # compute prior probabilities:
    start_prob = numpy.zeros((n_comps,))
    for i, u in enumerate(u_labels):
        start_prob[i] = numpy.count_nonzero(labels == u)
    # normalize prior probabilities
    start_prob = start_prob / start_prob.sum()

    # compute transition matrix:
    transmat = numpy.zeros((n_comps, n_comps))
    for i in range(labels.shape[0]-1):
        transmat[int(labels[i]), int(labels[i + 1])] += 1
    # normalize rows of transition matrix:
    for i in range(n_comps):
        transmat[i, :] /= transmat[i, :].sum()

    means = numpy.zeros((n_comps, n_feats))
    for i in range(n_comps):
        means[i, :] = numpy.matrix(features[:, 
                                   numpy.nonzero(labels == 
                                                 u_labels[i])[0]].mean(axis=1))

    cov = numpy.zeros((n_comps, n_feats))
    for i in range(n_comps):
        #cov[i,:,:] = numpy.cov(features[:,numpy.nonzero(labels==u_labels[i])[0]])  # use this lines if HMM using full gaussian distributions are to be used!
        cov[i, :] = numpy.std(features[:, numpy.nonzero(labels == 
                                                        u_labels[i])[0]], 
                              axis=1)

    return start_prob, transmat, means, cov