SSB - Sum of squares between clusters

Question

I got a little confused with the squares and the sums. As far as I know, the variance or total sum of squares (TSS) is smth like

$\sum_{i}^{n} (x_i - \bar x)^2$

and the sum of squares within (SSW) is

$\sum_{j}^{K} \sum_{i}^{n} (x_i - c_j)^2 \qquad i \in C_j$ where k ist the number of clusters

and that

$TSS = SSW + SSB$

Correct so far?

I therefore can do $TSS - SSW = SSB$

But what is the direct way to get $SSB$ from a given codebook? Preferablly I'd like to know, how to get there in numpy/scipy cause reading the equations is kind of hard for me..

import numpy as np
from scipy.cluster.vq import vq

# X.shape[0] -> observations
# X.shape[1] -> features
partition, euc_distance_to_centroids = vq(X, codebook)

TSS = np.sum((X-X.mean(0))**2)      
SSW = np.sum(euc_distance_to_centroids**2)
SSB = TSS - SSW

I think I'm missing the number of observations per cluster, when doing the SSB

>>> X
array([[ 2.,  4.,  2.],
       [ 1.,  3.,  1.],
       [ 3.,  4.,  2.],
       [ 2.,  3.,  2.],
       [ 1.,  5.,  5.]])
>>> codebook
array([[ 1.  ,  3.  ,  1.  ],
       [ 2.33,  3.67,  2.  ],
       [ 1.  ,  5.  ,  5.  ]])
>>> TSS
14.800000000000001
>>> SSW
1.3333333333333333
>>> SSB
13.466666666666667
>>> ((X.mean(0)-codebook)**2).sum() # How do I put the "num_clust_obs" in here?
12.542222222222223                  # Obviously not correct..

Answer 1

Could not find a shorter way, than this one

import numpy as np
from scipy.cluster.vq import vq

X = np.array([[ 2.,  4.,  2.],
              [ 1.,  3.,  1.],
              [ 3.,  4.,  2.],
              [ 2.,  3.,  2.],
              [ 1.,  5.,  5.]])

codebook = np.array([[ 1.  ,  3.  ,  1.  ],
                     [ 2.33,  3.67,  2.  ],
                     [ 1.  ,  5.  ,  5.  ]])

partition, euc_distance_to_centroids = vq(X, codebook)

TSS = np.sum((X-X.mean(0))**2)
SSW = np.sum(euc_distance_to_centroids**2)
SSB = TSS - SSW

# The 'direct' way
B = []
c = X.mean(0)
for i in range(partition.max()+1):
    ci = X[partition == i].mean(0)
    B.append(np.bincount(partition)[i]*np.sum((ci - c)**2))
SSB_ = np.sum(B)

print(TSS, SSW, SSB, SSB_)

prints

14.8 1.3334 13.4666 13.4666666667