This is a fix for the issue mentioned in #1708. It seems the -infs in the output were caused due to over-approximation for the bernoulli distribution. If I understand correctly, the num_gauss_hermite_points controls the level of approximation using the gaussian quadrature. If this number is high, some of the values during the computation of densities in variational_expectations is approximated as 1.0 (p). Hence, when you compute the log density for bernoulli here, it results in log(1-p) = log(1-1) = -inf.

I don't think this is a bug, but just a result of over-approximation as specified in the test. I re-run the test several times using 20 and it always passes. Please let me know if this seems reasonable or if you have any other thoughts/suggestions!

Thanks!