Gauss-Laguerre quadrature I am trying to compute this integral:
$$
\int_{0}^{\infty}\prod_{k = 1}^{d}\left(1 - \,\mathrm{e}^{-a_{k}\,t}\right)
\,\mathrm{e}^{-t}\,\mathrm{d}t,\quad
\mbox{where}\quad a_{k} > 0, \forall\ k.
$$
I can compute this just fine for small values $d$, e.g., less than $100$, using a numerical Gauss-Laguerre quadrature. I am having trouble computing this accurately when $d$ gets larger.  Any suggestions on how to solve this problem ?. 
 A: $\newcommand{\angles}[1]{\left\langle\,{#1}\,\right\rangle}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\half}{{1 \over 2}}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\iff}{\Longleftrightarrow}
 \newcommand{\imp}{\Longrightarrow}
 \newcommand{\Li}[1]{\,\mathrm{Li}_{#1}}
 \newcommand{\ol}[1]{\overline{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\ul}[1]{\underline{#1}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$

With $\ds{a_{k} > 0\,,\ \forall\ k \in \braces{1,2,\ldots,d}}$:

\begin{align}
&\color{#f00}{\int_{0}^{\infty}\prod_{k = 1}^{d}\pars{1 - \expo{-a_{k}\,t}}
\expo{-t}\,\dd t} =
\int_{0}^{\infty}\prod_{k = 1}^{d}
\pars{a_{k}t\int_{0}^{1}\expo{-a_{k}\,t\,x_{k}}\,\dd x_{k}}\expo{-t}\,\dd t
\\[3mm] = &\
\pars{\prod_{k = 1}^{d}a_{k}}
\int_{0}^{1}\cdots\int_{0}^{1}\int_{0}^{\infty}t^{d}
\exp\pars{-\bracks{\sum_{j = 1}^{d}a_{j}x_{j} + 1}t}\,\dd t
\,\prod_{k = 1}^{d}\dd x_{k}
\\[3mm] = &\
\Gamma\pars{d + 1}\pars{\prod_{k = 1}^{d}a_{k}}\,\mathrm{f}_{d}\pars{\vec{a}}
\tag{1}
\\[5mm] &\ \mbox{where}\
\,\mathrm{f}_{d}\pars{\vec{a}} \equiv
\int_{0}^{1}\cdots\int_{0}^{1}
{\prod_{k = 1}^{d}\,\dd x_{k} \over \pars{\vec{a}\cdot\vec{x} + 1}^{d + 1}}
\\[2mm] & \mbox{and}\quad
\left\lbrace\begin{array}{rcl}
\ds{\vec{a}} & \ds{\equiv} & \ds{\pars{a_{1},a_{2}\,\ldots,a_{d}}}
\\[2mm]
\ds{\vec{x}} &  \ds{\equiv} &  \ds{\pars{x_{1},x_{2}\,\ldots,x_{d}}}
\end{array}\right.
\end{align}

It's still a 'hard' one. However, when $\ds{d}$ is 'very large' the main contribution arises 'from' $\ds{\vec{a}\cdot\vec{x} \gtrsim 0}$ such that
\begin{align}
\,\mathrm{f}_{d}\pars{\vec{a}} & \sim
\int_{0}^{1}\cdots\int_{0}^{1}
\exp\pars{-\bracks{d + 1}\vec{a}\cdot\vec{x}}
\prod_{k = 1}^{d}\,\dd x_{k}
\\[5mm] & =
\prod_{k = 1}^{d}\int_{0}^{1}
\exp\pars{-\bracks{d + 1}a_{k}x_{k}}\,\dd x_{k} =
\prod_{k = 1}^{d}{1 - \expo{-\pars{d + 1}a_{k}} \over \pars{d + 1}a_{k}}
\\[5mm] & =
{1 \over \pars{d + 1}^{d}}
\prod_{k = 1}^{d}{1 - \expo{-\pars{d + 1}a_{k}} \over a_{k}}
\end{align}
and $\ds{\pars{~\mbox{see expression}\ \pars{1}~}}$
\begin{align}
&\color{#f00}{\int_{0}^{\infty}\prod_{k = 1}^{d}\pars{1 - \expo{-a_{k}\,t}}
\expo{-t}\,\dd t}
\\[3mm] \sim &\
\root{2\pi}d^{d + 1/2}\,\,\expo{-d}\pars{\prod_{k = 1}^{d}a_{k}}
{1 \over \pars{d + 1}^{d}}
\prod_{k = 1}^{d}{1 - \expo{-\pars{d + 1}a_{k}} \over a_{k}}\tag{2}
\\[3mm] & =
\root{2\pi}d^{1/2}{1 \over \pars{1 + 1/d}^{d}}\,\expo{-d} \sim
\color{#f00}{\root{2\pi}d^{1/2}\expo{-\pars{d + 1}}}
\end{align}

Note that
  $$
\bracks{1 - \expo{-\pars{d + 1}\min\braces{a_{k}}}\,\,}^{d} <
\prod_{k = 1}^{d}\bracks{1 - \expo{-\pars{d + 1}a_{k}}} <
\bracks{1 - \expo{-\pars{d + 1}\max\braces{a_{k}}}\,\,}^{d}
$$

A: Let we exploit $\int_{0}^{+\infty}e^{-kx}\,dx=\frac{1}{k}$ by expanding the product:
$$ I(\alpha_1,\ldots,\alpha_d)=\int_{0}^{+\infty}e^{-x}\prod_{k=1}^{d}\left(1-e^{-\alpha_k x}\right)\,dx \\= 1-\sum_{k=1}^{d}\frac{1}{\alpha_k+1}+\sum_{1\leq k_1<k_2\leq d}\frac{1}{\alpha_{k_1}+\alpha_{k_2}+1}-\ldots+(-1)^d\frac{1}{\alpha_1+\ldots+\alpha_d+1}\\=\color{blue}{\int_{0}^{1}\prod_{k=1}^{d}\left(1-x^{\alpha_k}\right)\,dx}$$
the last integral does not have a nice general closed form, but it is easy to approximate through Holder's inequality, since Euler's beta function gives:
$$ \int_{0}^{1}(1-x^{\alpha})^p\,dx = \frac{\Gamma(p+1)\,\Gamma\left(1+\frac{1}{a}\right)}{\Gamma\left(p+1+\frac{1}{a}\right)}.$$
A: The large number of factors will make the integrand look like a Gaussian near its peak. The maximum of the integrand is at $t = t^*$ satisfying the equation:
$$\sum_{j=1}^{d} \frac{a_j}{\exp(a_j t^*)-1} = 1$$
The second derivative of the logarithm of the integrand at the maximum is $-\sigma$ with
$$\sigma = \sum_{j=1}^{d} \frac{a_j^2\exp(a_j t^*)}{\left[\exp(a_j t^*)-1\right]^2}$$
The integrand near its peak is then approximately proportional to $\exp\left[-\frac{\sigma}{2}\left(t-t^*\right)^2\right]$. Using the values for $t^*$ and $\sigma$, you can set up an optimal Gauss-Hermite quadrature scheme. 
