A similar question was asked here and the given answer is perfect for a unidimensional integration.
I need to make bidimensional integration in R with a Gauss-Hermite quadrature: $$\int_{R^2} h(p1,p2) \phi(p1,p2) dp1 dp2$$ with $p1$ and $p2$ two parameters that follow a multivariate normal distribution of density $\phi(.)$ such as: $$p1,p2 \sim \mathcal{N_2} (\mu, \Sigma) $$ with $\mu = \left(\matrix{0\\ 1} \right)$ and $\Sigma = \left[\matrix{0.6&0.5\\0.5&0.25} \right]$.
I tried some packages in R but the results I got were very different from one another. I also tried the solution here but I'm not sure of the results I got.
I am looking for a simple solution, preferably based on gauss.quad.
Not being an expert, I have no idea of when to divide by $\pi$, if I should make a rescaling in my case and so on.
Thanks in advance for any help !
