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

def variational_expectations(self, Y, m, v, gh_points=None, Y_metadata=None):
        """
        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
        """

        if gh_points is None:
            gh_x, gh_w = self._gh_points()
        else:
            gh_x, gh_w = gh_points

        shape = m.shape
        m,v,Y = m.flatten(), v.flatten(), Y.flatten()

        #make a grid of points
        X = gh_x[None,:]*np.sqrt(2.*v[:,None]) + m[:,None]

        #evaluate the likelhood for the grid. First ax indexes the data (and mu, var) and the second indexes the grid.
        # broadcast needs to be handled carefully.
        logp = self.logpdf(X,Y[:,None], Y_metadata=Y_metadata)
        dlogp_dx = self.dlogpdf_df(X, Y[:,None], Y_metadata=Y_metadata)
        d2logp_dx2 = self.d2logpdf_df2(X, Y[:,None], Y_metadata=Y_metadata)

        #clipping for numerical stability
        #logp = np.clip(logp,-1e9,1e9)
        #dlogp_dx = np.clip(dlogp_dx,-1e9,1e9)
        #d2logp_dx2 = np.clip(d2logp_dx2,-1e9,1e9)

        #average over the gird to get derivatives of the Gaussian's parameters
        #division by pi comes from fact that for each quadrature we need to scale by 1/sqrt(pi)
        F = np.dot(logp, gh_w)/np.sqrt(np.pi)
        dF_dm = np.dot(dlogp_dx, gh_w)/np.sqrt(np.pi)
        dF_dv = np.dot(d2logp_dx2, gh_w)/np.sqrt(np.pi)
        dF_dv /= 2.

        if np.any(np.isnan(dF_dv)) or np.any(np.isinf(dF_dv)):
            stop
        if np.any(np.isnan(dF_dm)) or np.any(np.isinf(dF_dm)):
            stop

        if self.size:
            dF_dtheta = self.dlogpdf_dtheta(X, Y[:,None], Y_metadata=Y_metadata) # Ntheta x (orig size) x N_{quad_points}
            dF_dtheta = np.dot(dF_dtheta, gh_w)/np.sqrt(np.pi)
            dF_dtheta = dF_dtheta.reshape(self.size, shape[0], shape[1])
        else:
            dF_dtheta = None # Not yet implemented
        return F.reshape(*shape), dF_dm.reshape(*shape), dF_dv.reshape(*shape), dF_dtheta