code: def predictive_mean(self, mu, variance, Y_metadata=None):
          
          #conditional_mean: the edpected value of y given some f, under this likelihood
          fmin = -np.inf
          fmax = np.inf
          def int_mean(f,m,v):
              exponent = -(0.5/v)*np.square(f - m)
              #If exponent is under -30 then exp(exponent) will be very small, so don't exp it!)
              #If p is zero then conditional_mean will overflow
              assert v.all() > 0
              p = safe_exp(exponent)
  
              #If p is zero then conditional_variance will overflow
              if p < 1e-10:
                  return 0.
              else:
                  return self.conditional_mean(f)*p
          scaled_mean = [quad(int_mean, fmin, fmax,args=(mj,s2j))[0] for mj,s2j in zip(mu,variance)]
          mean = np.array(scaled_mean)[:,None] / np.sqrt(2*np.pi*(variance))
          return mean
docs: Quadrature calculation of the predictive mean: E(Y_star|Y) = E( E(Y_star|f_star, Y) )
  
          :param mu: mean of posterior
          :param sigma: standard deviation of posterior