Suppose I have a unimodal symmetric distribution of dimension $d$ > 3 such as the multivariate normal distribution ~ N(0, H), where H is a known $d$-dimensional covariance matrix.
I am interested in finding the Highest Density Region that covers 95% of the area under this distribution. Now, I understand the notion of Highest Density Regions in the univariate and bivariate case and how to proceed for obtaining these using for example the hdrcde package. However, I am interested in the following: How do I obtain the Highest Density Region given the distributional assumption of multivariate normality and a $d$-dimensional covariance matrix H where $d$=10 for example?
Idea how to proceed:
I am aware that in case of unimodal symmetric distributions such as the multivariate normal distribution considered in the question, the $\alpha$% HDR coincides with the $\alpha$% quantiles of the distribution. I could use a grid-based algorithm to approximate these quantiles for instance. The idea would be to define M grid points for each dimension such that relevant space is covered. This would would lead to M^$d$ joint grid values for $d$-dimensional distribution. Then I would evaluate the pdf of the multivariate normal density at each joint grid point. And then from this M^$d$ grid of densities I am not sure how to proceed in order to obtain $\alpha$-quantiles.
Any suggestions how to proceed are very welcome. Any other approach how to obtain High Density Regions for a $d$-dimensional normal distribution with $d>2$ are also more than welcome as I think M^$d$ function evaluations might become slow in higher dimensions. Thanking you in advance.