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.