I'm wondering if there is a good way to calculate the clustering criterion based on BIC formula, for a k-means output in R? I'm a bit confused as to how to calculate that BIC so that I can compare it with other clustering models. Currently I'm using the stats package implementation of k-means.

  • $\begingroup$ Note that this criterion is designed to be used with k-means. On clusterings obtained by different algorithms, it may be inappropriate (in particular for density based clustering algorithms) $\endgroup$ Apr 4 '13 at 18:56

To calculate the BIC for the kmeans results, I have tested the following methods:

  1. The following formula is from: [ref2] enter image description here

The r code for above formula is:

  k3 <- kmeans(mt,3)
  intra.mean <- mean(k3$within)
  k10 <- kmeans(mt,10)
  centers <- k10$centers
  BIC <- function(mt,cls,intra.mean,centers){
    x.centers <- apply(centers,2,function(y){
    sum1 <- sum(((mt-x.centers)/intra.mean)**2)
    sum1 + NCOL(mt)*length(unique(cls))*log(NROW(mt))

the problem is when i using the above r code, the calculated BIC was monotone increasing. what's the reason?

enter image description here

[ref2] Ramsey, S. A., et al. (2008). "Uncovering a macrophage transcriptional program by integrating evidence from motif scanning and expression dynamics." PLoS Comput Biol 4(3): e1000021.

  1. I have used the new formula from https://stackoverflow.com/questions/15839774/how-to-calculate-bic-for-k-means-clustering-in-r

    BIC2 <- function(fit){
    m = ncol(fit$centers)
        n = length(fit$cluster)
    k = nrow(fit$centers)
        D = fit$tot.withinss
    return(data.frame(AIC = D + 2*m*k,
                      BIC = D + log(n)*m*k))

This method given the lowest BIC value at cluster number 155. enter image description here

  1. using @ttnphns provided method, the corresponding R code as listed below. However, the problem is what the difference between Vc and V? And how to calculate the element-wise multiplication for two vectors with different length?

    BIC3 <- function(fit,mt){
    Nc <- as.matrix(as.numeric(table(fit$cluster)),nc=1)
    Vc <- apply(mt,2,function(x){
    V <- matrix(rep(apply(mt,2,function(x){
    LL = -Nc * colSums( log(Vc + V)/2 ) ##how to calculate this? elementa-wise multiplication for two vectors with different length?
    BIC = -2 * rowSums(LL) + 2*K*P * log(NRoW(mt))
  • 1
    $\begingroup$ Probably you were doing something different. It was stated in my "pseudocode" that Vc is P x K matrix and V was a column then propagated K times into the same sized matrix. So (point 4 in my answer) you can add Vc+V. Then take logarithm, divide by 2 and compute column sums. The resultant row vector multiply (value by value, i.e. elementwise) with row Nc. $\endgroup$
    – ttnphns
    Jun 30 '16 at 8:59
  • 1
    $\begingroup$ I've added formulas themselves in my answer, so you may compare and say if what you are doing is that or not that. $\endgroup$
    – ttnphns
    Mar 9 '19 at 12:18

I don't use R but here is a schedule which I hope will help you to compute the value of BIC or AIC clustering criteria for any given clustering solution.

This approach follows SPSS Algorithms Two-step cluster analysis (see the formulas there, starting from chapter "Number of clusters", then move to "Log-likelihood distance" where ksi, the log-likelihood, is defined). BIC (or AIC) is being computed based on the log-likelihood distance. I'm showing below computation for quantitative data only (the formula given in the SPSS document is more general and incorporates also categorical data; I'm discussing only its quantitative data "part"):

X is data matrix, N objects x P quantitative variables.
Y is column of length N designating cluster membership; clusters 1, 2,..., K.
1. Compute 1 x K row Nc showing number of objects in each cluster.
2. Compute P x K matrix Vc containing variances by clusters.
   Use denominator "n", not "n-1", to compute those, because there may be clusters with just one object.
3. Compute P x 1 column containing variances for the whole sample. Use "n-1" denominator.
   Then propagate the column to get P x K matrix V.
4. Compute log-likelihood LL, 1 x K row. LL = -Nc &* csum( ln(Vc + V)/2 ),
   where "&*" means usual, elementwise multiplication;
   "csum" means sum of elements within columns.
5. Compute BIC value. BIC = -2 * rsum(LL) + 2*K*P * ln(N),
   where "rsum" means sum of elements within row.
6. Also could compute AIC value. AIC = -2 * rsum(LL) + 4*K*P

Note: By default SPSS TwoStep cluster procedure standardizes all
quantitative variables, therefore V consists of just 1s, it is constant 1.
V serves simply as an insurance against ln(0) case.

AIC and BIC clustering criteria are used not only with K-means clustering. They may be useful for any clustering method which treats within-cluster density as within-cluster variance. Because AIC and BIC are to penalize for "excessive parameters", they unambiguously tend to prefer solutions with less clusters. "Less clusters more dissociated from one another" could be their motto.

There can be various versions of BIC/AIC clustering criteria. The one I showed here uses Vc, within-cluster variances, as the principal term of the log-likelihood. Some other version, perhaps better suited for k-means clustering, might base the log-likelihood on the within-cluster sums-of-squares.

The pdf version of the same SPSS document which I referred to.

And here is finally the formulae themselves, corresponding to the above pseudocode and the document; it is taken from the description of the function (macro) I've written for SPSS users. If you have any suggestions to improve the formulae please post a comment or an answer.

enter image description here

  • $\begingroup$ ttnphns, thanks for your response. I'm wondering if you could explain this with respect to the objective function which minimizes the within-cluster sum of squares? $\endgroup$ Apr 6 '13 at 20:00
  • $\begingroup$ You may see that the formula is almost the same as the second one in wiki. There, LL is based on error variance $\sigma_e^2$ which is Vc in my notation (the pooled within-cluster variance). Vc+V is used instead of Vc simply against the case of Vc=0 when a cluster has one object $\endgroup$
    – ttnphns
    Apr 6 '13 at 20:51
  • $\begingroup$ my head blows and i cant figure out how can i use this to stop my agglomerative hierarchical clustering . I am using it for document clustering record linkage problem $\endgroup$ Sep 4 '15 at 23:11
  • $\begingroup$ @Monster, There exist 100+ various internal clustering [validation] criterions. BIC is one of them. You do clustering to the end, saving cluster solutions, cluster membership variable on every step. Well, save only on last 10 or 20 steps because you probably don't want many small clusters. Compare the solutions by a clustering criterion and select 1-3 "best" ones. Compare them for interpretability. Done. See an example. $\endgroup$
    – ttnphns
    Sep 5 '15 at 10:06
  • $\begingroup$ @ttnphns the problem here i cant find any real data example of these so called 100+ internal clustering validation stuff. all i can find re some math equations which i can not understand. this also make me to disbelieve these algorithms exists in reality :D $\endgroup$ Sep 5 '15 at 11:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.