##### https://github.com/GPflow/GPflow
Tip revision: 6fd1a26
``````from __future__ import print_function, absolute_import
import itertools

import tensorflow as tf
import numpy as np

from gpflow import settings

def hermgauss(n):
x, w = np.polynomial.hermite.hermgauss(n)
x, w = x.astype(settings.np_float), w.astype(settings.np_float)
return x, w

def mvhermgauss(H, D):
"""
Return the evaluation locations 'xn', and weights 'wn' for a multivariate

The outputs can be used to approximate the following type of integral:
int exp(-x)*f(x) dx ~ sum_i w[i,:]*f(x[i,:])

:param H: Number of Gauss-Hermite evaluation points.
:param D: Number of input dimensions. Needs to be known at call-time.
:return: eval_locations 'x' (H**DxD), weights 'w' (H**D)
"""
gh_x, gh_w = hermgauss(H)
x = np.array(list(itertools.product(*(gh_x,) * D)))  # H**DxD
w = np.prod(np.array(list(itertools.product(*(gh_w,) * D))), 1)  # H**D
return x, w

def mvnquad(f, means, covs, H, Din, Dout=()):
"""
Computes N Gaussian expectation integrals of a single function 'f'
:param f: integrand function. Takes one input of shape ?xD.
:param means: NxD
:param covs: NxDxD
:param H: Number of Gauss-Hermite evaluation points.
:param Din: Number of input dimensions. Needs to be known at call-time.
:param Dout: Number of output dimensions. Defaults to (). Dout is assumed
to leave out the item index, i.e. f actually maps (?xD)->(?x*Dout).
"""
xn, wn = mvhermgauss(H, Din)
N = tf.shape(means)[0]

# transform points based on Gaussian parameters
chol_cov = tf.cholesky(covs)  # NxDxD
Xt = tf.matmul(chol_cov, tf.tile(xn[None, :, :], (N, 1, 1)), transpose_b=True)  # NxDxH**D
X = 2.0 ** 0.5 * Xt + tf.expand_dims(means, 2)  # NxDxH**D
Xr = tf.reshape(tf.transpose(X, [2, 0, 1]), (-1, Din))  # (H**D*N)xD