Calculating double integral with incomplete gamma function (numerically)

Question

I need to calculate a double integral on two variables ($b_0$ and $b_1$) in a function that includes a gamma incomplete function such as :

$\int_{\Bbb R} \frac{1}{\eta} (\lambda e^{b_{0} + b_{1}z_{1}})^{-1/\eta} \ \gamma\left(\frac{1}{\eta}, \lambda e^{b_{0} + b_{1}z_{1}} (t^*)^\eta \right) \frac{1} {(2\pi) \left| \boldsymbol{\Sigma}\right|^{1/2}}\;\; e^{ -\frac{1}{2}\left(\boldsymbol{b}-\boldsymbol{\mu}\right)^\top\boldsymbol{\Sigma}^{-1}\left(\boldsymbol{b}-\boldsymbol{\mu}\right) } db_0b_1$

The other variables in the function ($\lambda$, $\eta$...) are specified, $z$ is a vector. I am looking for a numerical resolution.

I already did some unsuccessful attempts in R.

Any advice to help me ? Thanks !

Answer 1

Two solutions (among ohers) may be used here:

• To use Gauss-Hermite quadrature since the second part of the expression is constitued of a bivariate normal density.
• The adaptIntegrate function from the cubaturepackage in R which is based on hypercubes