Gap Statistic Significance Testing with KMeans / KMedoids

Question

Relatively basic question here about FindClusters / ClusteringComponents. So I have noticed that Mathematica has an option for both KMeans / KMedoids as well as the Gap statistic for determining the optimal number of clusters. However, it seems that these can't be used together and that KMeans / KMedoids need the number of clusters entered beforehand. Both MATLAB and R allow these clustering methods to be used with the Gap statistic. Does anyone know of some Mathematica code floating around out there to implement this or is there already a way around this that I haven't spotted yet?

Thanks :)

Answer 1

I don't know of any pre-existing Mathematica code, but I had a look at the paper and gave it a go.

According to the paper, we are looking to compute $$ \rm{Gap}(k) = \rm{E}\{\log(W_k)\} - \log W_k $$

where $W_k$ is defined by $$ W_k = \sum_{r=1}^k\frac{1}{2n_r}D_r,\quad D_r=\sum_{ii'\in C_r}d_{ii'}. $$

Here, $C$ is a cluster computed by e.g. K-Means. $d$ is a distance metric. I believe there is a more efficient way of doing this with the distance to the mean of each cluster but I implemented $W_k$ the naive way:

w[x_] := Total[DistanceMatrix[x], 2]/(2 Length[x])
logwk[x_, k_] := Log@Total[w /@ FindClusters[x, k, Method -> "KMeans"]]

Actually, I implemented the log of $W_k$ since that is all we are ever going to need.

The second term in $\rm{Gap}(k)$ is simply logwk applied to the input data. The expectation, however, is actually the expectation of logwk applied to a reference distribution. A simple choice for the reference distribution is the uniform distribution. Basically, we're going to compare logwk applied to our data and logwk applied to the uniform distribution, if our data is clustered then for the right number of clusters $k$, $W_k$ will be smaller for our data than for the uniform distribution, because $W_k$ is a measurement of inter-cluster variance and is optimized by finding well-defined clusters.

Here is code to generate data from the reference distribution:

ref[x_] := Module[{xrange, yrange},
  {xrange, yrange} = CoordinateBounds[x];
  Transpose[{
    RandomReal[xrange, Length[x]],
    RandomReal[yrange, Length[x]]
    }]
  ]

We will estimate the expectation by taking the mean value of $W_k$ over $B$ realizations of ref[x]. Furthermore, because we are using a finite sample, we will introduce a correction using the standard deviation. All in all, the gap statistic looks like this:

gapStatistic[x_, k_, B_: 10.] := Module[{refValues, sd, mu},
  refValues = Table[logwk[ref[x], k], {B}];
  sd = StandardDeviation[refValues];
  mu = Mean[refValues];
  {mu - logwk[x, k], mu - logwk[x, k] - sd}
  ]

I return the value with and without the correction because both are needed in the final step, which is to check what the smallest $k$ is so that, $$ \rm{Gap}(k) > Gap(k+1) - s_{k+1}, $$ where $s$ is the correction.

gapTest[x_, kmax_, B_: 10.] := Module[{gaps, comps},
  gaps = Table[gapStatistic[x, k, B], {k, kmax}];
  LengthWhile[Rest@Partition[gaps, 2, 1], #[[1, 1]] < #[[2, 2]] &] + 2
  ]

I skipped the case with one cluster because it was giving me trouble (we will actually see this in the following example.) Other than that, it seems to be working. Let's do an example:

sampleData[k_] := Flatten[Table[RandomVariate[
    MultinormalDistribution[
     RandomReal[{-10, 10}, 2],
     IdentityMatrix[2]
     ], 50], {k}], 1]

SeedRandom[101];
data = sampleData[5];

ListPlot[data]

We would expect gapTest to say that 5 is the optimal number of clusters, and it does:

gapTest[data, 10]

5

Let's have a look at the curves for the left and right-hand side of inequality:

gaps = Table[gapStatistic[data, k, 10.], {k, 10}];
ListLinePlot[{
  gaps[[;; -2, 1]],
  gaps[[2 ;;, 2]]
  }]

We see that except for at $k=1$ it behaves exactly as we would expect, with the left-hand side becoming larger than the right-hand side at $k=5$.

The function is rather slow but I have no time to optimize it now, unfortunately.