code: def moments_match_ep(self,obs,tau,v,Y_metadata_i=None):
          
          #Compute first integral for zeroth moment.
          #NOTE constant np.sqrt(2*pi/tau) added at the end of the function
          mu = v/tau
          sigma2 = 1./tau
          #Lets do these for now based on the same idea as Gaussian quadrature
          # i.e. multiply anything by close to zero, and its zero.
          f_min = mu - 20*np.sqrt(sigma2)
          f_max = mu + 20*np.sqrt(sigma2)
  
          def int_1(f):
              return self.pdf(f, obs, Y_metadata=Y_metadata_i)*np.exp(-0.5*tau*np.square(mu-f))
          z_scaled, accuracy = quad(int_1, f_min, f_max)
  
          #Compute second integral for first moment
          def int_2(f):
              return f*self.pdf(f, obs, Y_metadata=Y_metadata_i)*np.exp(-0.5*tau*np.square(mu-f))
          mean, accuracy = quad(int_2, f_min, f_max)
          mean /= z_scaled
  
          #Compute integral for variance
          def int_3(f):
              return (f**2)*self.pdf(f, obs, Y_metadata=Y_metadata_i)*np.exp(-0.5*tau*np.square(mu-f))
          Ef2, accuracy = quad(int_3, f_min, f_max)
          Ef2 /= z_scaled
          variance = Ef2 - mean**2
  
          #Add constant to the zeroth moment
          #NOTE: this constant is not needed in the other moments because it cancells out.
          z = z_scaled/np.sqrt(2*np.pi/tau)
  
          return z, mean, variance
docs: Calculation of moments using quadrature
  
          :param obs: observed output
          :param tau: cavity distribution 1st natural parameter (precision)
          :param v: cavity distribution 2nd natural paramenter (mu*precision)