`groupings` is described on the doc string as a `list`, but I believe it is an `array`. 

The calls each time to `groupings == i` will get expensive if both `groupings` and `number_groups` can get large. A different numpythonic approach here would be:

```python
for idx, group_ix in enumerate(groupings == arange(number_groups)[:, newaxis]):
  diffs = distances[group_ix]
  s_W += (diffs ** 2).sum() / unique_n[idx]
```

note, the motivation also to use `array.sum()` instead of `sum` is that python's builtin `sum` is much slower.