Starting with a remark. I have a $P$-dimensional integral I'm evaluating for which I know the "right" answer and have developed an implementation that yields back the right answer (based on unit testing with other implementations and simulation). I'm boiling down many pages of work simply to try and ask my question and I hope I have done so coherently, but will add/edit this to reflect more information if needed.
The integral is of the following general form
$$ \widehat{\boldsymbol{\theta}} = \frac{\int \dots \int \boldsymbol{\theta} L (\boldsymbol{k}|\boldsymbol{\theta}) f (\boldsymbol{\theta}\boldsymbol{|}\boldsymbol{\mu},\boldsymbol{\Sigma})d\boldsymbol{\theta}}{\int \dots \int L(\boldsymbol{k}|\boldsymbol{\theta})f(\boldsymbol{\theta}\boldsymbol{|}\boldsymbol{\mu},\boldsymbol{\Sigma})d\boldsymbol{\theta}} $$
where $L(\boldsymbol{k|\theta})$ is a likelihood function consisting of a product of binomial densitites and $f(\boldsymbol{\theta})$ is multivariate gaussian with mean vector $\boldsymbol{\mu}$ and covariance $\boldsymbol{\Sigma}$. In this case the vector $\boldsymbol{\mu} = 0$ for all elements.
My implementation evaluates the integral via summation over an $Q^P \times P$ grid of nodes and weights where $Q$ is the number of quadrature points chosen in one dimension. Let me now define $\theta_q$ and $w_{q}$ as the nodes and weights, respectively, from the one-dimensional normalized gauss hermite rule (e.g., https://juliastats.org/MixedModels.jl/v1.1/GaussHermite/#Evaluating-the-weights-and-abscissae-1). My grid of nodes is based on the cartesian product of these one-dimensional (original) nodes and the weights ($w^*_{q}$) are the product of weights over all dimensions for the corresponding nodes in each row of this grid.
My question now is whether my implementation is an accurate reflection of the Liu and Pierce work (https://www.jstor.org/stable/2337136) which states we can rewrite this integral to include using a ratio of normal densities to better sample in the region needed, or adapt into the region.
My estimator is $$ \widehat{\theta}_p = \frac{\sum^Q_{q=1}\theta_{p,q} L(\boldsymbol{k}_1|\theta_{1,q}),\ldots, L(\boldsymbol{k}_p|\theta_{P,q})f(\boldsymbol{\theta}\boldsymbol{|}\boldsymbol{\mu},\boldsymbol{\Sigma})w^*_q}{\sum^Q_{q=1}L(\boldsymbol{k}_1|\theta_{1,q}), \ldots, L(\boldsymbol{k}_p|\theta_{P,q})f(\boldsymbol{\theta}\boldsymbol{|}\boldsymbol{\mu},\boldsymbol{I})w^*_q}\tag{1} $$
Here, the density in the numerator is over the vector of points $\boldsymbol{\theta} = (\theta_{1,q}, \ldots,\theta_{p,q})'$ and the same in the denominator but here I am using $\boldsymbol{I}$ (identity matrix), which through experimentation and trial is yielding back good results.
Some implementations I have read rotate the nodes via the cholesky factor of the covariance matrix, $\boldsymbol{\Sigma} = \boldsymbol{L}\boldsymbol{L}'$ as in $\boldsymbol{L}\boldsymbol{\theta}$ but I am not doing so.
My implementation in (1) yields back the right answer (again from comparison in to other implementations) and I want to verify my application of the Liu and Pierce theory is accurate.
Thank you in advance for advice, clarifications, and corrections.