code: def psicomputations(variance, Z, variational_posterior, return_psi2_n=False):
      
      # here are the "statistics" for psi0, psi1 and psi2
      # Produced intermediate results:
      # psi0    N
      # psi1    NxM
      # psi2    MxM
      mu = variational_posterior.mean
      S = variational_posterior.variance
      gamma = variational_posterior.binary_prob
  
      psi0 = (gamma*(np.square(mu)+S)*variance).sum(axis=-1)
      psi1 = np.inner(variance*gamma*mu,Z)
      psi2 = np.inner(np.square(variance)*(gamma*((1-gamma)*np.square(mu)+S)).sum(axis=0)*Z,Z)+tdot(psi1.T)
  
      return psi0, psi1, psi2
docs: Compute psi-statistics for ss-linear kernel