I am implementing Welford's algorithm to compute mean and variance on a series of images.

Broadcasting does not work.

The algorithm for scalars courtesy of Wikipedia

def online_variance(data):
    n = 0
    mean = 0.0
    M2 = 0.0

    for x in data:
        n += 1
        delta = x - mean
        mean += delta/n
        delta2 = x - mean
        M2 += delta*delta2

    if n < 2:
        return float('nan')
    else:
        return M2 / (n - 1)

For numpy arrays I am forced to replace mean += delta/n by mean = mean + delta/n

The error I get is non-broadcastable output operand with shape () doesn't match the broadcast shape (256,256,3). My images have shape 256x256x3 (RGB).