code: def update_statistics(Y, P, beta=0.9):
    
  
    (n, d) = Y.shape
    
  
    A = np.expand_dims(Y, axis=1) * np.expand_dims(Y, axis=0)
    assert(A.shape == (n, n, d))
    
    Q = np.mean(A, axis=2)
    Q[np.where(Q == 0.)] = 0.000001
    assert(P.shape == Q.shape)
  
    return beta*P + (1-beta)*Q
docs: Args
    ----
    Y: 
        2d array whose columns encode the 
        activity of the output units
    P:
        2d array encoding the pairwise average 
        activity of the output units  
  
    Returns
    -------
        The updated average activities