def predictive_variance(self, mu,variance, predictive_mean=None, Y_metadata=None):
<EOL>normalizer = np.sqrt(<NUM_LIT:2>*np.pi*variance)<EOL>fmin_v = -np.inf<EOL>fmin_m = np.inf<EOL>fmin = -np.inf<EOL>fmax = np.inf<EOL>from ..util.misc import safe_exp<EOL>def int_var(f,m,v):<EOL><INDENT>exponent = -(<NUM_LIT:0.5>/v)*np.square(f - m)<EOL>p = safe_exp(exponent)<EOL>if p < <NUM_LIT>:<EOL><INDENT>return <NUM_LIT:0.><EOL><DEDENT>else:<EOL><INDENT>return self.conditional_variance(f)*p<EOL><DEDENT><DEDENT>scaled_exp_variance = [quad(int_var, fmin_v, fmax,args=(mj,s2j))[<NUM_LIT:0>] for mj,s2j in zip(mu,variance)]<EOL>exp_var = np.array(scaled_exp_variance)[:,None] / normalizer<EOL>if predictive_mean is None:<EOL><INDENT>predictive_mean = self.predictive_mean(mu,variance)<EOL><DEDENT>predictive_mean_sq = predictive_mean**<NUM_LIT:2><EOL>def int_pred_mean_sq(f,m,v,predictive_mean_sq):<EOL><INDENT>exponent = -(<NUM_LIT:0.5>/v)*np.square(f - m)<EOL>p = np.exp(exponent)<EOL>if p < <NUM_LIT>:<EOL><INDENT>return <NUM_LIT:0.><EOL><DEDENT>else:<EOL><INDENT>return self.conditional_mean(f)**<NUM_LIT:2>*p<EOL><DEDENT><DEDENT>scaled_exp_exp2 = [quad(int_pred_mean_sq, fmin_m, fmax,args=(mj,s2j,pm2j))[<NUM_LIT:0>] for mj,s2j,pm2j in zip(mu,variance,predictive_mean_sq)]<EOL>exp_exp2 = np.array(scaled_exp_exp2)[:,None] / normalizer<EOL>var_exp = exp_exp2 - predictive_mean_sq<EOL>return exp_var + var_exp<EOL>
Approximation to the predictive variance: V(Y_star)

The following variance decomposition is used:
V(Y_star) = E( V(Y_star|f_star)**2 ) + V( E(Y_star|f_star) )**2

:param mu: mean of posterior
:param sigma: standard deviation of posterior
:predictive_mean: output's predictive mean, if None _predictive_mean function will be called.