MATLAB Answers

Need fast quadrature method for use within root-finding routine

2 views (last 30 days)
Bob Johnson
Bob Johnson on 9 Aug 2011
Edited: John D'Errico on 9 Oct 2020
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,

  1 Comment

Mike Hosea
Mike Hosea on 10 Aug 2011
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?

Sign in to comment.

Answers (1)

John D'Errico
John D'Errico on 9 Oct 2020
Edited: John D'Errico on 9 Oct 2020
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.


Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by