Need fast quadrature method for use within root-finding routine

Hi I need to find the root of a function that contains an expectation over a three-dimensional multivariate normal distribution, where the dimensions are all independent. This algorithm will be embedded within an estimation routine over a parameter space. So, the quadrature has to be fast. The matlab routines are too unreliable to embed within an estimation routine, and in some cases too slow. I usually use the multivariate normal quadrature routine in Comp Econ. But, there is a discrete change in behavior in of the dimension. I know the location of the discrete point, so I could split the integral up in that dimension, but does anyone have an idea about a good, fast quadrature routine over a truncated normal? The number of nodes needed for good approximation with Newton Coates Rules makes those methods insufficient.

Thanks for your help. All the best,

Bob

[Mike Hosea]

QUADGK and QUAD2D are too unreliable or too slow? With QUAD2D it's sometimes important not to integrate over discontinuities and areas of non-smoothness (because in some cases it attracts a very fine mesh, so it's often better to split the integral into pieces), but apart from that, I would have expected them to do fairly well. Could you provide an example that illustrates the problem with these codes?

Answer 1

[John D'Errico]

Way too late for this answer to be useful. In fact, way too late even for it to have been useful for me, when I had to solve the same problem, back in maybe 1978 or so?

But what I would do now, is tessellate the region of interest, dissecting into tetrahedrae. In the case where the mode of the Gaussian is inside the region, I would insure that point was one of the nodes of said dissection. If the Gaussian mode is not inside the region, then just dissect it arbitrarily. You can dissect a cube into either 5 or 6 tetrahedrae - your choice. The dissection into 6 pieces is nice in a way, in that they are all similar tetrahedrae.

Then inside every simplex, use a Gaussian style quadrature, designed to work on a simplex. I recall finding some that apply to 3-d. Perhaps better yet would be to use a simple rule, and then split each simplex into a set of smaller similar simplexes. Use a Richardson estrapolation to measure convergence and speed it up, subdividing only those pieces as needed.

I would add the above scheme works reasonably well for 2-d problems.