def variational_expectations(self, Y, m, v, gh_points=None, Y_metadata=None):
if gh_points is None:<EOL><INDENT>gh_x, gh_w = self._gh_points()<EOL><DEDENT>else:<EOL><INDENT>gh_x, gh_w = gh_points<EOL><DEDENT>shape = m.shape<EOL>m,v,Y = m.flatten(), v.flatten(), Y.flatten()<EOL>X = gh_x[None,:]*np.sqrt(<NUM_LIT>*v[:,None]) + m[:,None]<EOL>logp = self.logpdf(X,Y[:,None], Y_metadata=Y_metadata)<EOL>dlogp_dx = self.dlogpdf_df(X, Y[:,None], Y_metadata=Y_metadata)<EOL>d2logp_dx2 = self.d2logpdf_df2(X, Y[:,None], Y_metadata=Y_metadata)<EOL>F = np.dot(logp, gh_w)/np.sqrt(np.pi)<EOL>dF_dm = np.dot(dlogp_dx, gh_w)/np.sqrt(np.pi)<EOL>dF_dv = np.dot(d2logp_dx2, gh_w)/np.sqrt(np.pi)<EOL>dF_dv /= <NUM_LIT><EOL>if np.any(np.isnan(dF_dv)) or np.any(np.isinf(dF_dv)):<EOL><INDENT>stop<EOL><DEDENT>if np.any(np.isnan(dF_dm)) or np.any(np.isinf(dF_dm)):<EOL><INDENT>stop<EOL><DEDENT>if self.size:<EOL><INDENT>dF_dtheta = self.dlogpdf_dtheta(X, Y[:,None], Y_metadata=Y_metadata) <EOL>dF_dtheta = np.dot(dF_dtheta, gh_w)/np.sqrt(np.pi)<EOL>dF_dtheta = dF_dtheta.reshape(self.size, shape[<NUM_LIT:0>], shape[<NUM_LIT:1>])<EOL><DEDENT>else:<EOL><INDENT>dF_dtheta = None <EOL><DEDENT>return F.reshape(*shape), dF_dm.reshape(*shape), dF_dv.reshape(*shape), dF_dtheta<EOL>
Use Gauss-Hermite Quadrature to compute

   E_p(f) [ log p(y|f) ]
   d/dm E_p(f) [ log p(y|f) ]
   d/dv E_p(f) [ log p(y|f) ]

where p(f) is a Gaussian with mean m and variance v. The shapes of Y, m and v should match.

if no gh_points are passed, we construct them using defualt options