def METHOD_NAME(self, D, means, labels, K, R):
    """Computes the Bayesian Information Criterion."""
    D = vq.whiten(D)
    Rn = D.shape[0]
    M = D.shape[1]
    if R == K:
        return 1
    # Maximum likelihood estimate (MLE)
    mle_var = 0
    for k in range(len(means)):
        X = D[np.argwhere(labels == k)]
        X = X.reshape((X.shape[0], X.shape[-1]))
        for x in X:
            mle_var += distance.euclidean(x, means[k])
            #print x, means[k], mle_var
    mle_var /= float(R - K)
    # Log-likelihood of the data
    l_D = - Rn/2. * np.log(2*np.pi) - (Rn * M)/2. * np.log(mle_var) - \
        (Rn - K) / 2. + Rn * np.log(Rn) - Rn * np.log(R)
    # Params of BIC
    p = (K-1) + M * K + mle_var
    #print "BIC:", l_D, p, R, K
    # Return the bic
    return l_D - p / 2. * np.log(R)