def predictive_mean(self, mu, variance, Y_metadata=None):
<EOL>fmin = -np.inf<EOL>fmax = np.inf<EOL>def int_mean(f,m,v):<EOL><INDENT>exponent = -(<NUM_LIT:0.5>/v)*np.square(f - m)<EOL>assert v.all() > <NUM_LIT:0><EOL>p = safe_exp(exponent)<EOL>if p < <NUM_LIT>:<EOL><INDENT>return <NUM_LIT:0.><EOL><DEDENT>else:<EOL><INDENT>return self.conditional_mean(f)*p<EOL><DEDENT><DEDENT>scaled_mean = [quad(int_mean, fmin, fmax,args=(mj,s2j))[<NUM_LIT:0>] for mj,s2j in zip(mu,variance)]<EOL>mean = np.array(scaled_mean)[:,None] / np.sqrt(<NUM_LIT:2>*np.pi*(variance))<EOL>return mean<EOL>
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