Multivariate integral involving a Gaussian and a product of t-Student densities

Question

Let $N \ge 1$ be an integer and let $A$ be a positive definite $N\times N$ matrix and let $b$ be a $N\times 1$ real vector.Let $\nu > 2$ and define $n:= (\nu+1)/2$.

We are interested in deriving the distribution of a sample covariance matrix where the observations are being drawn from a specific non-Gaussian distribution. Here we decide to replace the Gaussian by a t-Student distribution. This would be an analogue of the Wishart distribution where the population is t-Student. Now, it turns out that computing the Laplace transform of such a matrix leads to the following multivariate integral below.

\begin{eqnarray} {\mathfrak f}_\nu^{(A,b)} &:=& \int\limits_{{\mathbb R}^N} \frac{e^{-\frac{1}{2} z^{T} \cdot A \cdot z + b^T \cdot z} }{\sqrt{(2\pi)^N \cdot \det(A^{-1})}} \cdot \left( \prod\limits_{i=1}^N \frac{d z_i}{\left(1 + \frac{z_i^2}{\nu}\right)^n} \right) \tag{1} \end{eqnarray}

The integral above is a multi-variate generalization of a certain uni-variate integral as given in the other post. We define a new matrix which is produced from the original matrix by multiplying the diagonal elements by inverses of integration variables, i.e ${\tilde A}(\eta) := \left( A_{i,j}\cdot (1 + \delta_{i,j} (\frac{1}{\eta_i}-1) ) \right)_{i,j=1}^N$.

Now, by using the representation of a t-Student density as a continuous mixture of normal distributions with their inverse variances being Gamma distributed and then swapping the order of integration , doing the multi-variate Gaussian integral and simplifying the result we came up with the following result below.

\begin{eqnarray} &&{\mathfrak f}_\nu^{(A,b)} = e^{\frac{\nu}{2} \cdot Tr[A]} \cdot \sqrt{\det[A]} \cdot \left(\prod\limits_{i=1}^N (\frac{\nu}{2} A_{i,i})^n \right) \cdot \frac{1}{((n-1)!)^N} \cdot \int\limits_{(0,1)^N} \left[\prod\limits_{i=1}^N \frac{1}{\eta_i^2} \left( \frac{1-\eta_i}{\eta_i} \right)^{n-1} \cdot e^{-\frac{\nu}{2} \frac{A_{i,i}}{\eta_i} } \right] % \frac{\exp\left( \frac{1}{2} b^T \cdot [{\tilde A}(\eta)]^{-1} \cdot b\right)}{\sqrt{\det {\tilde A}(\eta) }} % \left[\prod\limits_{i=1}^N d \eta_i \right] \tag{2} \end{eqnarray}

The Mathematica code snippet below confirms that the result is correct. Here we go:

(*Dimension equals NN.*)
NN = 3; Clear[z, A, b, n, \[Eta]];
A = RandomReal[{-1, 1}, {NN, NN}]; A = 
 Transpose[A] . DiagonalMatrix[RandomReal[{1/2, 1}, NN]] . 
  A; A += +IdentityMatrix[NN];
PositiveDefiniteMatrixQ[A]
b = RandomReal[{1, 2}, {NN, 1}];
\[Nu] = RandomReal[{3, 5}];
n = (\[Nu] + 1)/2;

1/Sqrt[(2 \[Pi])^NN Det[A]^-1]
  NIntegrate[
  Exp[-1/2  Sum[A[[i, j]] z[i]  z[j], {i, 1, NN}, {j, 1, NN}] + 
     Sum[b[[i]] z[i], {i, 1, NN}]]  Product[1/(1 + z[i]^2/\[Nu])^
    n, {i, 1, NN}], 
  Evaluate[
   Sequence @@ Table[{z[i], -Infinity, +Infinity}, {i, 1, NN}]]]

Exp[\[Nu]/2  Tr[A]]  Sqrt[Det[A]]  Product[(\[Nu]/2  A[[i, i]])^
  n, {i, 1, NN}] 1/((n - 1)!)^NN  NIntegrate[
  Product[1/\[Eta][i]^2 ((1 - \[Eta][i])/\[Eta][i])^(n - 1)
      Exp[-\[Nu]/2 A[[i, i]]/\[Eta][i]], {i, 1, NN}] Exp[
    1/2  Flatten[
       b] . (Inverse[
         Table[If[i == j, 1/\[Eta][i], 1] A[[i, j]], {i, 1, NN}, {j, 
           1, NN}]] . Flatten[b])]/Sqrt[
   Det[Table[
     If[i == j, 1/\[Eta][i], 1] A[[i, j]], {i, 1, NN}, {j, 1, NN}]]], 
  Evaluate[Sequence @@ Table[{\[Eta][i], 0, 1}, {i, 1, NN}]]]

$(2)$." />

---

Having said all that above my question is how do we use the stationary phase approximation to come up with an approximate closed form for our quantity as $\nu$ is big?

Another question is how do we evaluate that integral in a different way?

Answer 1

We essentially proceed in the same way as in the $N=1$ case. Firstly we rewrite the first term in square brackets in the right hand side in $(2)$. Here you go:

\begin{eqnarray} &&{\mathfrak f}_\nu^{(A,b)} = e^{\frac{\nu}{2} \cdot Tr[A]} \cdot \sqrt{\det[A]} \cdot \left(\prod\limits_{i=1}^N (\frac{\nu}{2} A_{i,i})^n \right) \cdot \frac{1}{((n-1)!)^N} \cdot\\ && \int\limits_{(0,1)^N} % \underline{ \left[ \prod\limits_{i=1}^N e^{-n \frac{A_{i,i}}{\eta_i} + n \log(\frac{1-\eta_i}{\eta_i})} \right] } \cdot \underbrace{ \left[ \prod\limits_{i=1}^N \frac{e^{\frac{1}{2} \frac{A_{i,i}}{\eta_i} }}{\eta_i(1-\eta_i)} \right] % % \frac{\exp\left( \frac{1}{2} b^T \cdot [{\tilde A}(\eta)]^{-1} \cdot b\right)}{\sqrt{\det {\tilde A}(\eta) }} } % \left[\prod\limits_{i=1}^N d \eta_i \right] \tag{3} \end{eqnarray}

Now, the phase function is the argument of the exponential in the underlined term. It is easy to see that the stationary point of that phase function reads $ \eta_* = \left( \frac{A_{i,i}}{1 + A_{i,i}} \right)_{i=1}^N $.

Now, we do two things. Firstly, we expand the argument of the exponential in the underlined term about the stationary point and then we expand the remaining part of the integrand about that point. In both cases we use expansion to the second order. We define ${\tilde {\tilde A}} := A + 1_N + \sum\limits_{\xi=1}^\infty \frac{(-1)^\xi}{\sqrt{n}^\xi} \cdot D_\xi$ with $D_\xi:= diag \left( \frac{\eta_i^\xi}{(1+A_{i,i})^{\xi-1}}\right)_{i=1}^N$ and $B:=(A+1_N)^{-1}$ and after some fairly simple manipulations we obtain the following result below.

\begin{eqnarray} &&{\mathfrak f}_\nu^{(A,b)} \simeq % % \left(\underline{ \sqrt{2\pi} (n-\frac{1}{2})^n \frac{e^{-(n-\frac{1}{2})}}{(n-1)! \sqrt{n} } }\right)^N \cdot % \sqrt{\det(A)} \cdot \frac{\exp\left( \frac{1}{2} b^T \cdot [A+1_N]^{-1} \cdot b\right)}{\sqrt{\det(A+1_N)}} \cdot \\ && % \int\limits_{-(1+A_{1,1}) \sqrt{n}}^{\frac{1+A_{1,1}}{A_{1,1}} \sqrt{n}} \cdots \int\limits_{-(1+A_{N,N}) \sqrt{n}}^{\frac{1+A_{N,N}}{A_{N,N}} \sqrt{n}} % \left( \prod\limits_{i=1}^N \frac{e^{-\frac{1}{2} \eta_i^2}}{\sqrt{2\pi}} \right) \cdot \\ && \left( \prod\limits_{i=1}^N e^{ % -\frac{\eta _i^3 \left(A_{i,i}-2\right)}{3 \sqrt{n} \left(A_{i,i}+1\right)} % -\frac{\eta _i^4 \left(A_{i,i}^2-2 A_{i,i}+3\right)}{4 n \left(A_{i,i}+1\right){}^2} % +O(\frac{1}{n^{3/2}}) } \right) \cdot \\ && \left( \prod\limits_{i=1}^N [1+ \frac{-3+A_{i,i}}{2(1+A_{i,i})} \cdot \frac{\eta_i}{\sqrt{n}} + \frac{17 - 2 A_{i,i} + 5 A_{i,i}^2}{8(1+A_{i,i})^2} \cdot (\frac{\eta_i}{\sqrt{n}})^2 + O(\frac{1}{n^{3/2}}) ] \right) \cdot \\ && \left( \right. \\ && \left. 1+ \right. \\ && \left. \frac{1}{2 \sqrt{n}} \cdot (\text{Tr}\left[B.D_1\right]+b^T.B.D_1.B.b)+ \right. \\ &&\left. % \frac{1}{8 n} \cdot (\left(\text{Tr}\left[B.D_1\right]+b^T.B.D_1.B.b\right){}^2-4 \left(\text{Tr}\left[B.D_2\right]+b^T.B.D_2.B.b-b^T.B.D_1.B.D_1.B.b\right)+2 \text{Tr}\left[B.D_1.B.D_1\right]) \right. \\ &&\left. + O(\frac{1}{n^{3/2}}) \right. \\ &&\left. \right)\\ && % % % \left[\prod\limits_{i=1}^N d \eta_i \right] \tag{4} \end{eqnarray}

Now, as always, the code below verifies the result numerically for $N=3$ and for randomly chose $A,b$ and $\nu$ (subject to $\nu \ge 10$). Here we go:

In[2461]:= (*Dimension equals NN.*)
NN = 3; Clear[z, A, b, n, \[Eta]];
A = RandomReal[{-1, 1}, {NN, NN}]; A = 
 Transpose[A] . DiagonalMatrix[RandomReal[{1/2, 1}, NN]] . 
  A; A += +IdentityMatrix[NN];
PositiveDefiniteMatrixQ[A]
b = RandomReal[{1, 2}, {NN, 1}];
\[Nu] = RandomReal[{10, 15}];
n = (\[Nu] + 1)/2;

1/Sqrt[(2 \[Pi])^NN Det[A]^-1]
  NIntegrate[
  Exp[-1/2 Sum[A[[i, j]] z[i] z[j], {i, 1, NN}, {j, 1, NN}] + 
     Sum[b[[i]] z[i], {i, 1, NN}]] Product[1/(1 + z[i]^2/\[Nu])^
    n, {i, 1, NN}], 
  Evaluate[
   Sequence @@ Table[{z[i], -Infinity, +Infinity}, {i, 1, NN}]]]

Exp[\[Nu]/2 Tr[A]]  Sqrt[Det[A]]  Product[(\[Nu]/2 A[[i, i]])^
  n, {i, 1, NN}] 1/((n - 1)!)^NN  NIntegrate[
  Product[Exp[-n A[[i, i]]/\[Eta][i] + 
      n  Log[(1 - \[Eta][i])/\[Eta][i]]], {i, 1, NN}] Product[
    Exp[(1/2) A[[i, i]]/\[Eta][i]]/(\[Eta][i] (1 - \[Eta][i])), {i, 1,
      NN}] Exp[
    1/2 Flatten[
       b] . (Inverse[
         Table[If[i == j, 1/\[Eta][i], 1] A[[i, j]], {i, 1, NN}, {j, 
           1, NN}]] . Flatten[b])]/Sqrt[
   Det[Table[
     If[i == j, 1/\[Eta][i], 1] A[[i, j]], {i, 1, NN}, {j, 1, NN}]]], 
  Evaluate[Sequence @@ Table[{\[Eta][i], 0, 1}, {i, 1, NN}]]]
(*Here we applied the stationary phase approximation.*)

Exp[\[Nu]/2  Tr[A]]  Sqrt[Det[A]]  Product[(\[Nu]/2 A[[i, i]])^
  n, {i, 1, NN}] 1/((n - 1)!)^NN  NIntegrate[
  Product[Exp[
     n  Log[1/A[[i, i]]] - n  (1 + A[[i, i]]) - (
      n (1 + A[[i, i]])^4 (-(A[[i, i]]/(1 + A[[i, i]])) + \[Eta][
          i])^2)/(2 A[[i, i]]^2)], {i, 1, NN}] Product[(
     E^(1/2 (1 + A[[i, i]])) (1 + A[[i, i]])^2)/A[[i, i]] + (
     E^(1/2 + 
       1/2 A[[i, i]]) (-3 + A[[i, i]]) (1 + 
        A[[i, i]])^3 (-(A[[i, i]]/(1 + A[[i, i]])) + \[Eta][i]))/(
     2 A[[i, i]]^2) + (
     E^(1/2 + 
       1/2 A[[i, i]]) (1 + A[[i, i]])^4 (17 - 2 A[[i, i]] + 
        5 A[[i, i]]^2) (-(A[[i, i]]/(1 + A[[i, i]])) + \[Eta][i])^2)/(
     8 A[[i, i]]^3), {i, 1, 
     NN}] (Exp[
      1/2 Flatten[
         b] . (Inverse[
           Table[If[
              i == j, (1 + A[[i, i]])/
               A[[i, i]] - ((1 + 
                  A[[i, i]])^2 (-(A[[i, i]]/(1 + A[[i, i]])) + \[Eta][
                   i]))/A[[i, 
                i]]^2 + ((1 + 
             
                    A[[i, i]])^3 (-(A[[i, i]]/(
                   1 + A[[i, i]])) + \[Eta][i])^2)/A[[i, i]]^3, 1] A[[
              i, j]], {i, 1, NN}, {j, 1, NN}]] . 
          Flatten[b])]/(\[Sqrt]Det[
        Table[If[
           i == j, (1 + A[[i, i]])/
            A[[i, i]] - ((1 + 
               A[[i, i]])^2 (-(A[[i, i]]/(1 + A[[i, i]])) + \[Eta][
                i]))/A[[i, 
             i]]^2 + ((1 + 
               A[[i, i]])^3 (-(A[[i, i]]/(1 + A[[i, i]])) + \[Eta][
                i])^2)/A[[i, i]]^3, 1] A[[i, j]], {i, 1, NN}, {j, 1, 
          NN}]])), 
  Evaluate[Sequence @@ Table[{\[Eta][i], 0, 1}, {i, 1, NN}]]]


 Sqrt[Det[A]]  (Sqrt[2 \[Pi]] (n - 1/2 )^
   n Exp[-(n - 1/2)]/((n - 1)! Sqrt[n]))^NN  NIntegrate[
  Product[Exp[-((\[Eta][i])^2/2 )]/Sqrt[2 \[Pi]], {i, 1, NN}] Product[
    1 +  (-3 + A[[i, i]])  /(
      2  (1 + A[[i, i]])^1) (\[Eta][i]/Sqrt[
       n]) +   (17 - 2 A[[i, i]] + 5  A[[i, i]]^2) /(
      8  (1 + A[[i, i]])^2) (\[Eta][i]/Sqrt[n])^2, {i, 1, 
     NN}] (Exp[
      1/2  Flatten[
         b] . (Inverse[
           Table[If[i == j, 
             1 + A[[i, i]] - (\[Eta][i]/Sqrt[n]) + (\[Eta][i]/Sqrt[
               n])^2/ (1 + A[[i, i]]), A[[i, j]]], {i, 1, NN}, {j, 1, 
             NN}]] . Flatten[b])]/(\[Sqrt]Det[
        Table[If[i == j, 
          1 + A[[i, i]] - (\[Eta][i]/Sqrt[n]) + (\[Eta][i]/Sqrt[
            n])^2/ (1 + A[[i, i]]), A[[i, j]]], {i, 1, NN}, {j, 1, 
          NN}]])), 
  Evaluate[
   Sequence @@ 
    Table[{\[Eta][i], - (1 + A[[i, i]])^1 Sqrt[n], (1 + A[[i, i]])^1/
       A[[i, i]] Sqrt[n]}, {i, 1, NN}]]]



Out[2463]= True

Out[2467]= {1.38337}

During evaluation of In[2461]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

Out[2468]= 1.38337

Out[2469]= 1.54514

Out[2470]= 1.54514

Overall conclusion:

It is not hard to see that when expanding the integrand in $(4)$ in powers of $1/\sqrt{n} $ we will obtain a multivariate expansion in the vector $\eta$ with firstly linear terms, then quadratic and possibly higher order terms. The zeroth order term factorizes into a product of Gaussian cumulative functions (CDFs) whereas the higher order terms will be always expressed in terms of derivatives of those CDFs. This basically concludes the whole thing.

Answer 2

Here we provide a different way of computing the integral in question. This is based on the idea that the (un-normalized) Student-t density $\rho_\nu(z) := \left(1+ z^2/\nu\right)^{-\nu/2-1/2} $ can be expanded in a Taylor series about $\nu=0$. We will be truncating that expansion at some threshold ${\mathfrak n} \ge 1$. To be specific let us define some auxiliary quantities first. Here we go:

\begin{eqnarray} P(z) &:=& \left( \left. \frac{1}{n!} \frac{d^n}{d x^n} (1+z^2 x)^{-\frac{1}{2}(1+\frac{1}{x})} \right|_{x=0} e^{\frac{z^2}{2}} \right)_{n=0}^{{\mathfrak n}} \\ &=& \left( \begin{array}{c} 1 \\ \frac{1}{4} z^2 \left(z^2-2\right) \\ \frac{1}{96} z^4 \left(3 z^4-28 z^2+36\right) \\ \frac{1}{384} z^6 \left(z^6-22 z^4+116 z^2-120\right) \\ \vdots \end{array} \right) \tag{1} \end{eqnarray}

The quantities above are the coefficients at inverse powers of $\nu$ in the expansion of $\rho_\nu(z)/\rho_\infty(z)$.

Now we define other quantities below. Here we go:

\begin{eqnarray} Q_\theta(\vec{z}) &:=& \sum\limits_{\begin{array}{lll} \xi_1+\cdots + \xi_N = \theta \\ \xi_1 \ge 0, \cdots, \xi_N \ge 0 \end{array} } \prod\limits_{i=1}^N P_{\xi_i}(z_i) \\ &=& \left( \begin{array}{c} 1 \\ \sum\limits_{i=1}^N P_1(z_i) \\ \sum\limits_{i=1}^N P_2(z_i) + \sum\limits_{1 \le i < j \le N} P_1(z_i) P_1(z_j) \\ \sum\limits_{i=1}^N P_3(z_i) + \sum\limits_{1 \le i < j \le N} [ P_2(z_i) P_1(z_j) + P_1(z_i) P_2(z_j)] + \sum\limits_{1 \le i < j < k \le N} P_1(z_i) P_2(z_j) P_3(z_k) \\ \vdots \end{array} \right) \tag{2} \end{eqnarray}

for $\theta = 0,\cdots, {\mathfrak n}$.

The quantities above are the coefficients at inverse powers of $\nu$ in the expansion of $\prod_{i=1}^N \rho_\nu(z_i)$.

With those definitions at hand the result reads:

\begin{eqnarray} {\mathfrak f}^{(A,b)}_\nu = \frac{1}{\sqrt{\det(1_N + A^{-1})}} \cdot \sum\limits_{\theta=0}^{\mathfrak n} \nu^{-\theta} Q_\theta(\frac{\vec{\partial}}{\partial \vec{b}}) \cdot e^{\frac{1}{2} \vec{b}^T \cdot (A+1_N)^{-1} \cdot \vec{b} } \tag{3} \end{eqnarray}

---

We have tested formula $(3)$ numerically .We took $N=3$ and $10 \le \nu \le 15$ and ${\mathfrak n} = 6$ and we conclude that the result is accurate to a five decimal digits precision. Here we go:

(*Dimension equals NN. Method II.*)
NN = 3; Clear[z, x, A, b, n, z, zz, Q, B];
nmax = 6; mprec = 50;
P[z_] = Simplify@
   CoefficientList[
    Normal[Series[(1 + z^2  x)^(-1/2  1/x - 1/2), {x, 0, nmax}]] E^(
     z^2/2), x];

A = RandomReal[{-1, 1}, {NN, NN}, WorkingPrecision -> mprec]; A = 
 Transpose[A] . 
  DiagonalMatrix[
   RandomReal[{1/2, 1}, NN, WorkingPrecision -> mprec]] . 
  A; A += +IdentityMatrix[NN];
PositiveDefiniteMatrixQ[A]
b = RandomReal[{1, 2}, {NN}, WorkingPrecision -> mprec];
\[Nu] = RandomReal[{10, 15}, WorkingPrecision -> mprec];
n = (\[Nu] + 1)/2; \[Epsilon] = 1/Sqrt[n];

B = (A + IdentityMatrix[NN]);
BI = Inverse[(A + IdentityMatrix[NN])];
zz = Table[z[i], {i, 1, NN}];
Q[zz_List] = 
  Table[Sum[
    Product[P[
       z[i]][[1 + 
        If[i < NN, \[Xi][i], \[Theta] - 
          Sum[\[Xi][j], {j, 1, NN - 1}]]]] , {i, 1, NN}], 
    Evaluate[
     Sequence @@ 
      Table[{\[Xi][i], 
        0, \[Theta] - Sum[\[Xi][j], {j, 1, i - 1}]}, {i, 1, 
        NN - 1}]]], {\[Theta], 0, nmax}];

(*Calculations.*)

1/Sqrt[(2 \[Pi])^NN Det[A]^-1]
  NIntegrate[
  Exp[-1/2  Sum[A[[i, j]] z[i]  z[j], {i, 1, NN}, {j, 1, NN}] + 
     Sum[b[[i]] z[i], {i, 1, NN}]]  Product[1/(1 + z[i]^2/\[Nu])^
    n, {i, 1, NN}], 
  Evaluate[
   Sequence @@ Table[{z[i], -Infinity, +Infinity}, {i, 1, NN}]]]

Clear[z, \[Xi]];
1/Sqrt[(2 \[Pi])^NN Det[A]^-1]
  NIntegrate[
  Exp[-1/2  Sum[B[[i, j]] z[i]  z[j], {i, 1, NN}, {j, 1, NN}] + 
     Sum[b[[i]] z[i], {i, 1, NN}]]  Sum[
    Q[zz][[1 + \[Theta]]] \[Nu]^-\[Theta], {\[Theta], 0, nmax}], 
  Evaluate[
   Sequence @@ Table[{z[i], -Infinity, +Infinity}, {i, 1, NN}]]]

Clear[z];
Clear[\[Xi], \[Xi]\[Xi]];
SS = Expand[
   Table[Q[zz][[1 + \[Theta]]] \[Nu]^-\[Theta], {\[Theta], 0, 
     nmax}]];
\[Xi]\[Xi] = Table[\[Xi][j], {j, 1, NN}];
ll = Table[{(# /. z[i_] :> 1), Exponent[#, zz]} & /@ (List @@ 
      SS[[1 + \[Theta]]]), {\[Theta], 0, nmax}];
ll[[1]] = {{1, ConstantArray[0, NN]}};

expTrms = ConstantArray[0, {nmax + 1}];
Do[
  t0 = TimeUsed[];
  expTrms[[1 + \[Theta]]] = 
   1/Sqrt[Det[A]^-1] 1/Sqrt[Det[B]]
      Nest[Sum[
        ll[[1 + \[Theta], i, 1]]  D[#, 
          Evaluate[
           Sequence @@ 
            Table[{\[Xi][j], ll[[1 + \[Theta], i, 2, j]]}, {j, 1, 
              NN}]]], {i, 1, Length[ll[[1 + \[Theta]]]]}] &, 
      Exp[1/2 \[Xi]\[Xi] . (BI . \[Xi]\[Xi])], 1] /. \[Xi][j_] :> 
     b[[j]];
  Print["\[Theta] = ", \[Theta], " done in time=", TimeUsed[] - t0];
  , {\[Theta], 0, nmax}];
MatrixForm[#] & /@ {expTrms, Accumulate[expTrms]}