Use
self.variances = np.sum(gamma * (img[:, None] - self.means) ** 2, axis=0) / N_k
instead of
self.variances = np.sum(gamma * (img - self.means) ** 2) / N_k
import numpy as np
N, K = 10, 20
gamma = np.random.random((N, K))
means = np.random.random(K)
N_k = np.random.random(K)
img = np.random.random(N)
expected = np.array([ np.sum([gamma[n][k] * (img[n] - means[k]) ** 2 for n in range(0, N)]) / N_k[k] for k in range(0, K) ])
result = np.sum(gamma * (img[:, None] - means) ** 2, axis=0) / N_k
assert np.allclose(result, expected)
Note that img - self.means subtracts corresponding values elementwise.
img has shape (N,) and self.means has shape (K,), so subtracting one from the other is possible if N == K, but raises an ValueError if N != K since then the shapes are incompatible. Since you didn't get an error, N must be equal to K.
Since you want to compute a 2D array of values, first make img a 2D array of shape (N, 1) by adding a new axis: img[:, None].
Then take advantage of
broadcasting so
that when you subtract the 1D array self.means, it will broadcast to a 2D
array of compatible shape (N, K). Broadcasting by default adds new axes on
the left. Thus the array of shape (K,) can automatically broadcast to an
array of shape (N, K). (This also explains why we needed to use img[:, None] to explicitly add a new axis on the right.)
Now (img[:, None] - self.means) will be a 2D array of shape (N, K).
Also note that it is important to specify axis=0 when calling np.sum because the summation is done over the first axis of length N. Since Python uses 0-based indexing the first axis corresponds to axis=0.
If you don't specify an axis, then np.sum sums over all axes by default.
