When programming machine learning programs that contain integrals, it is often pretty hard to both do the integrals themselves and also to translate them into the computer. It typically involves “vectorising” the calculations, which is often not straightforward. Because of this, it’s hard to get it right the first time.
If the integrals are like the ones that appear in Bayesian inference, however, it is very easy to write a Monte Carlo estimator for them:
$$ \int f(x) p(x) dx = \sum_{n=1}^N f(x_i) \text{ for } x_i \sim p(x_i) $$
Ideally then, one would write a unit test that runs your integral code and compares it to a slower and less precise version, the Monte Carlo estimator. But, precisely because of its Monte Carlo nature, you cannot just compare the output to the output from the integral code.