1
$\begingroup$

As Wikipedia teaches us https://en.wikipedia.org/wiki/Owen%27s_T_function the Owen's T function $T(h,a)$ defines a probability of a bivariate event $X>h$ and $0<Y<a X$ where $X,Y$ are standard, independent Gaussian random variables.

Now in the context of question Multivariate gaussian integral over positive reals a necessity appeared to deal with a slightly more general quantity. \begin{equation} T(h,a,b):= {\bf P}\left(X>h \quad \wedge \quad a X+b > Y > 0 \left. \right| X = N(0,1) , Y=N(0,1) \right) \end{equation} We have shown that : \begin{eqnarray} &&T(h,a,b)= \int\limits_h^\infty \frac{\exp(-1/2 \xi^2)}{\sqrt{2\pi}} \frac{1}{2} Erf(\frac{a \xi+b}{\sqrt{2}}) d\xi \quad (i1)\\ &&= \int\limits_0^a \frac{e^{-\frac{b^2}{2}-b h \xi -\frac{1}{2} h^2 \left(\xi ^2+1\right)}}{2 \pi \left(\xi ^2+1\right)} d\xi - \frac{b}{2\sqrt{2}\sqrt{\pi}} \int\limits_0^a \frac{\xi e^{-\frac{b^2}{2 \xi ^2+2}} \text{erfc}\left(\frac{\xi (b+h \xi )+h}{\sqrt{2} \sqrt{\xi ^2+1}}\right)}{\left(\xi ^2+1\right)^{3/2}} d\xi + \frac{1}{4} \text{erf}\left(\frac{b}{\sqrt{2}}\right) \text{erfc}\left(\frac{h}{\sqrt{2}}\right) \quad (i2) \end{eqnarray}

{a, b, h} = RandomReal[{0, 1}, 3, WorkingPrecision -> 50]; b = 0;
NIntegrate[
 Exp[-x^2/2]/Sqrt[2 Pi] 1/2 Erf[(a x + b)/Sqrt[2]], {x, h, Infinity}, 
 WorkingPrecision -> 20]
NIntegrate[(E^(-(b^2/2) - xi b h - 1/2 (1 + xi^2) h^2)) /(
   2 (1 + xi^2) \[Pi]) - 
   b  /(2 Sqrt[2] Sqrt[ \[Pi]]) (
    xi  Erfc[(h + xi (b + xi h))/(Sqrt[2] Sqrt[1 + xi^2])])/ ((1 + 
      xi^2)^(3/2)) E^(-(b^2/(2 + 2 xi^2))), {xi, 0, a}, 
  WorkingPrecision -> 20] + Erfc[h/Sqrt[2]] Erf[b/Sqrt[2]] 1/4

enter image description here

Update: Let $A_j \in {\mathbb R}$ for $j=1,\cdots,3$ and let $x\in {\mathbb R}$. Then we have: \begin{eqnarray} T(A_1 x, A_2, A_3 x) = \frac{1}{2\pi} \left(\arctan(A_2)-\arctan(A_2+\frac{A_3}{A_1})-\arctan(\frac{A_1+A_2 A_3+A_2^2 A_1}{A_3})\right) + \frac{1}{4} erf[\frac{A_3 x}{\sqrt{2} \sqrt{1+A_2^2}}] + T(A_1 x, \frac{A_2 A_1+A_3}{A_1})+ T(\frac{A_3 x}{\sqrt{1+A_2^2}},\frac{A_1+A_2 A_3 + A_2^2 A_1}{A_3}) \quad (ii) \end{eqnarray}

This identity comes from differentiating both sides with respect to $x$ then using the definition of the generalized Owen's T function to evaluate the derivative on the right hand side and having done this integrating both sides with respect to $x$ again.

Let us present the proof of that in detail. Firstly we define $f(x) := T[A_1 x, A_2, A_3 x]$. Now we compute the derivative using the chain rule. We have: \begin{eqnarray} \frac{d }{d x} f(x) &=& \partial_1 T[A_1 x, A_2 , A_3 x] \cdot A_1 + \partial_3 T[A_1 x, A_2, A_3 x] \cdot A_3 \\ &=& - \left. \rho(h) \frac{1}{2} erf[\frac{a h+b}{\sqrt{2}}] \right|_{\begin{array}{r} h=A_1 x \\ a=A_2 \\ b=A_3 x \end{array}}\cdot A_1 + \left.\frac{1}{\sqrt{1+a^2}} \frac{1}{2} erf[\frac{h+a^2 h+a b}{\sqrt{2} \sqrt{1+a^2}}] \rho(\frac{b}{1+a^2})\right|_{\begin{array}{r} h=A_1 x \\ a=A_2 \\ b=A_3 x \end{array}} \cdot A_3 \\ &=& -\rho(A_1 x) \frac{1}{2} erf[\frac{A_1 A_2 + A_3}{\sqrt{2}} x] \cdot A_1 + \frac{1}{\sqrt{1+A_2^2}} \rho(\frac{A_3 x}{\sqrt{1+A_2^2}}) \frac{1}{2} erfc[\frac{A_1+A_2 A_3 +A_1 A_2^2}{\sqrt{2} \sqrt{1+A_2^2}} x] \cdot A_3 \end{eqnarray}

Now we integrate. We have: \begin{eqnarray} f(x)- f(0) &=& - \int\limits_0^x \rho(A_1 \xi) \frac{1}{2} erf[\frac{A_1 A_2 + A_3}{\sqrt{2}} \xi] d\xi \cdot A_1 + \\ &&\frac{1}{\sqrt{1+A_2^2}}\int\limits_0^x \rho(\frac{A_3 \xi}{\sqrt{1+A_2^2}}) \frac{1}{2} erfc[\frac{A_1+A_2 A_3 +A_1 A_2^2}{\sqrt{2} \sqrt{1+A_2^2}} \xi] d\xi \cdot A_3 \\ f(x) - \frac{1}{2\pi} \arctan(A_2) &=& - \frac{1}{2\pi} \arctan\left( \frac{A_1 A_2+A_3}{A_1}\right) + T(A_1 x, \frac{A_1 A_2 + A_3}{A_1}) + \\ && \frac{1}{4} erf\left( \frac{A_3}{\sqrt{2} \sqrt{1+A_2^2}} x\right) +\\ &&-\frac{1}{2\pi} \arctan\left( \frac{A_1+A_2 A_3 + A_1 A_2^2}{A_3}\right) + T\left( \frac{A_3}{\sqrt{1+A_2^2}} x, \frac{A_1+A_2 A_3 + A_1 A_2^2}{A_3}\right) \end{eqnarray} where in the second line we used the results from An integral involving error functions and a Gaussian and the definition of the Owen's T function. This completes the proof.

(*A certain derivative. Used in Q869502.nb*)
T[h_, a_, b_] := 
  NIntegrate[(E^(-(b^2/2) - xi b h - 1/2 (1 + xi^2) h^2)) /(
     2 (1 + xi^2) \[Pi]) - 
     b  /(2 Sqrt[2] Sqrt[ \[Pi]]) (
      xi  Erfc[(h + xi (b + xi h))/(Sqrt[2] Sqrt[1 + xi^2])])/ ((1 + 
        xi^2)^(3/2)) E^(-(b^2/(2 + 2 xi^2))), {xi, 0, a}, 
    WorkingPrecision -> 20] + Erfc[h/Sqrt[2]] Erf[b/Sqrt[2]] 1/4;
{A1, A2, A3} = RandomReal[{-1, 1}, 3, WorkingPrecision -> 50];
u = Range[0, 1, 1/100];
mT = Interpolation[Transpose[{u, T[A1 u, A2, A3 u]}]];
u =.; u = RandomReal[{0, 1}, WorkingPrecision -> 50];
mT'[u]
-rho[A1 u] 1/2 Erf[(A1 A2 + A3)/Sqrt[2] u] A1 + 
 1/Sqrt[1 + A2^2]
   rho[(A3 u)/Sqrt[1 + A2^2]] 1/
   2 Erfc[(A1 + A2 A3 + A1 A2^2)/(Sqrt[2] Sqrt[1 + A2^2]) u] A3

T[A1 u, A2, A3 u]
1/(2 Pi) (ArcTan[A2] - ArcTan[(A2 A1 + A3)/A1] - 
    ArcTan[(A1 + A2 A3 + A2^2 A1)/A3]) + 
 1/4 Erf[(A3 u)/(Sqrt[2] Sqrt[1 + A2^2])] + 
 OwenT[A1 u, (A2 A1 + A3)/A1] + 
 OwenT[A3/Sqrt[1 + A2^2] u, (A1 + A2 A3 + A2^2 A1)/A3]
1/(2 Pi) (-ArcTan[A3/((A1 + A2 A3 + A2^2 A1))] - 
    ArcTan[(A1 + A2 A3 + A2^2 A1)/A3]) + 
 1/4 Erf[(A3 u)/(Sqrt[2] Sqrt[1 + A2^2])] + 
 OwenT[A1 u, (A2 A1 + A3)/A1] + 
 OwenT[A3/Sqrt[1 + A2^2] u, (A1 + A2 A3 + A2^2 A1)/A3]
-1/(2 Pi) Pi/2 (Sign[A3/((A1 + A2 A3 + A2^2 A1))]) + 
 1/4 Erf[(A3 u)/(Sqrt[2] Sqrt[1 + A2^2])] + 
 OwenT[A1 u, (A2 A1 + A3)/A1] + 
 OwenT[A3/Sqrt[1 + A2^2] u, (A1 + A2 A3 + A2^2 A1)/A3]


-(1/4) Sign[A3/((A1 + A2 A3 + A2^2 A1))] + 
 1/4 Erf[(A3 u)/(Sqrt[2] Sqrt[1 + A2^2])] + 
 OwenT[A1 u, (A2 A1 + A3)/A1] + 
 OwenT[A3/Sqrt[1 + A2^2] u, (A1 + A2 A3 + A2^2 A1)/A3]

enter image description here

Now by both taking $x=1$ and replacing $A_1$,$A_2$ and $A_3$ by $h$, $a$ and $b$ in $(ii)$ we express the generalized Owen's T function through Owen's T function itself. We have: \begin{eqnarray} T(h,a,b) = \frac{1}{2\pi} \left(\arctan(a)-\arctan(a+\frac{b}{h})-\arctan(\frac{h+a b+a^2 h}{b})\right) + \frac{1}{4} erf[\frac{b}{\sqrt{2(1+a^2)}}] + T\left( h,\frac{a h+b}{h}\right) + T\left( \frac{b}{\sqrt{1+a^2}},\frac{h+a b+a^2 h}{b}\right) \end{eqnarray}

As a sanity check we look at the limit $b$ going to zero. We have: \begin{eqnarray} \lim_{b \rightarrow 0_+} T(h,a,b) &=& \frac{1}{2\pi} \left(\arctan(a)-\arctan(a)-\frac{\pi}{2} sign(h))\right) + 0 + T(h,a) + \frac{1}{4} sign(h) \\ &=& T(h,a) \end{eqnarray} as it should be.

As another sanity check we look at the case $a=\imath$. Going back to the calculations of the derivative above we have: \begin{eqnarray} \frac{d}{d x} f(x)= -\phi(A_1 x) \frac{1}{2} erf(\frac{A_1 A_2+A_3}{\sqrt{2}} x) A_1 + \frac{1}{2\pi \imath x} \exp(-\frac{1}{2} x^2 (2 A_1 \imath A_3+A_3^2)) \end{eqnarray} where we used the asymptotic expansion for the complementary error function given in https://en.wikipedia.org/wiki/Error_function#Complementary_error_function . Now we take a number $M$ such that $1< M$ and we integrate the above from unity to $M$ and we get: \begin{eqnarray} f(1)-f(M)= \left.\left( T(A_1 \cdot \xi,A_2+\frac{A_3}{A_1}) + \frac{1}{4\pi \imath} Ei(-\frac{1}{2}(1+2\imath \frac{A_1}{A_3})(\xi A_3)^2\right)\right|_{\xi=M}^{\xi=1} \end{eqnarray} where $Ei()$ is the exponential integral. Now it turns out that as $M\rightarrow \infty$ both $f(M)$ and $T(\dots M,\dots)$ tend to zero and \begin{equation} \lim\limits_{M\rightarrow \infty} \frac{1}{4 \pi \imath} Ei((a+\imath b) M) = sign(b) \cdot \frac{1}{4} \cdot 1_{a<0} + \infty \cdot 1_{a>0} \end{equation} Defining $b:=b_1+\imath b_2$ and taking $h>0$ this gives the final result: \begin{eqnarray} &&T(h,\imath , b) = \\ &&\left\{ \begin{array}{rr} T(h,\imath + \frac{b}{h}) + \frac{1}{4\pi \imath} Ei(\frac{1}{2}(-b_1^2+b_2^2+2 b_2 h-2\imath b_1(b_2+h))) + sign(b_1(b_2+h)) \cdot \frac{1}{4} & \mbox{if $b_2<0$ and $-b_1^2 + b_2^2+2 b_2 h <0$} \\ \infty & \mbox{otherwise} \end{array} \right. \end{eqnarray}

My question is the following. Has this quantity been ever analyzed in the literature before?

$\endgroup$
0
$\begingroup$

This is not an answer to the question above but instead it is a generalization of the results above. Define $\vec{a}:=(a_j)_{j=1}^d \in {\mathbb R}_+^d$ and let us define as $T^{(d)}(h,\vec{a})$ the probability of a following multivariate event $X>h$ and $0< Y_j < a_j X$ for $j=1,\cdots,d$ where $X$ and $\left( Y_j \right)_{j=1}^d$ are standard, independent Gaussian random variables.

Now take another vector $\vec{b}:=(b_j)_{j=1}^d \in {\mathbb R}_+^d$ and define a slightly more general quantity: \begin{eqnarray} T^{(d)}(h,\vec{a},\vec{b}) &:=& P\left( X > h \quad \wedge \quad \begin{array}{rrr} a_1 X + b_1 >& Y_1& >0 \\ \vdots &\vdots& \vdots \\ a_d X + b_d >& Y_d& >0 \end{array}\left. \right| \begin{array}{rrr} X&=&N(0,1)\\Y_1 &=&N(0,1)\\ &\vdots& \\ Y_d&=&N(0,1) \end{array} \right)\\ &=& \int\limits_h^\infty \rho(\xi) \left[\prod\limits_{i=1}^d \frac{1}{2} erf(\frac{a_i \xi+b_i}{\sqrt{2}})\right] d\xi \end{eqnarray} In what follows we will prove that if $d\le 2$ then the quantity $T^{(d)}(h,\vec{a},\vec{b})$ reduces to elementary functions and to $T^{(d)}(h,\vec{a})$ only.

As in the question above we consider a following quantity $T^{(d)}(h \cdot x, \vec{a}, \vec{b} \cdot x)$ which we differentiate with respect to $x$. We have: \begin{eqnarray} &&\frac{d}{d x} T^{(d)}(h \cdot x, \vec{a}, \vec{b} \cdot x)=\\ &&-h \cdot \rho(h \cdot x) \prod\limits_{i=1}^d \frac{1}{2} erf(\frac{a_i h x+b_i x}{\sqrt{2}})+\\ &&\sum\limits_{i=1}^d \frac{b_i}{\sqrt{2 \pi}} \int\limits_{h \cdot x}^\infty e^{-\frac{1}{2} (a_i \xi + b_i x)^2}\rho(\xi) \left[\prod\limits_{j=1,j\neq i}^d \frac{1}{2} erf(\frac{a_j h x+b_j x}{\sqrt{2}})\right] d\xi \end{eqnarray} What we do now is to simplify the second term on the right hand side, i.e. absorb the exponential into the Gaussian density and extract a constant prefactor. After that we integrate the identity above over $x$ from zero to unity. The result reads: \begin{eqnarray} &&T^{(d)}(h,\vec{a},\vec{b}) - T^{(d)}(0,\vec{a},\vec{0})=\\ &&-T^{(d)}(0,\vec{a}+\frac{1}{h} \vec{b},\vec{0}) + T^{(d)}(h,\vec{a}+\frac{1}{h} \vec{b},\vec{0}) + \\ &&\sum\limits_{i=1}^d \int\limits_0^{\frac{b_i}{\sqrt{1+a_i^2}}} \rho(x) \cdot T^{(d-1)}\left([a_i + \frac{h}{b_i}(1+a_i^2)]x,\frac{(a_j)_{j=1,j\neq i}^d}{\sqrt{1+a_i^2}},\frac{(b_j(1+a_i^2)-b_i a_i a_j)_{j=1,j\neq i}^d}{b_i \sqrt{1+a_i^2}} \right) dx \end{eqnarray} this clearly gave us a recurrence relation for the quantity in question subject to $T^{(0)}(h,\vec{a},\vec{b})=T^{(0)}(h)= 1/2 erfc(h/\sqrt{2})$.

Now we state the result for $d=2$. Firstly we define auxiliary quantities: \begin{eqnarray} \delta&:=& h^2+(a_1 h+b_1)^2+(a_2 h+b_2)^2\\ \delta_1&:=&h(1+a_1^2+a_2^2) + a_1 b_1 + a_2 b_2\\ \delta_2&:=&1+a_1^2+a_2^2\\ \hline\\ (m_1,m_2)&:=& (b_1(1+a_2^2)-a_1 a_2 b_2,b_2(1+a_1^2)-a_1 a_2 b_1)\\ (n_1,n_2)&:=&(h+h a_1^2+a_1 b_1,h+h a_2^2+a_2 b_2)\\ (o_1,o_2)&:=&(h a_1+b_1,h a_2 +b_2)\\ (p_1,p_2)&:=&(1+a_1^2,1+a_2^2)\\ (k_1,k_2)&:=&(\frac{\sqrt{p_1} \delta_1}{m_2},\frac{\sqrt{p_2} \delta_1}{m_1})\\ (l_1,l_2)&:=&(\frac{m_1}{\sqrt{p_2 \delta_2}},\frac{m_2}{\sqrt{p_1 \delta_2}}) \end{eqnarray} Then the result reads: \begin{eqnarray} &&4 \pi T^{(2)}(h,\vec{a},\vec{b})=\\ &&\arctan(\frac{a_1 a_2}{\sqrt{\delta_2}}) - \arctan(\frac{o_1 o_2}{h \sqrt{\delta}})+\\ &&\arctan(\frac{m_2}{\sqrt{\delta_2} b_1}) + \arctan(\frac{m_1}{\sqrt{\delta_2} b_2})+\\ &&\arctan(\frac{b_1 o_2}{n_1 \sqrt{\delta}}) + \arctan(\frac{b_2 o_1}{n_2 \sqrt{\delta}})+\\ &&\arctan(\frac{b_1 \delta_1}{\sqrt{m_2^2 \delta}}) + \arctan(\frac{b_2 \delta_1}{\sqrt{m_1^2 \delta}})+\\ &&\left( \arctan(\frac{a_2}{\sqrt{p_1}}) - \arctan(\frac{\sqrt{p_1} o_2}{n_1})-\arctan(k_1)\right) \cdot erf(\frac{b_1}{\sqrt{2 p_1}})+\\ &&\left( \arctan(\frac{a_1}{\sqrt{p_2}}) - \arctan(\frac{\sqrt{p_2} o_1}{n_2})-\arctan(k_2)\right) \cdot erf(\frac{b_2}{\sqrt{2 p_2}})+\\ &&2\pi \left( T(\frac{n_1}{\sqrt{p_1}},\frac{\sqrt{p_1} o_2}{n_1}) + T(l_2,k_1)\right)\cdot erf(\frac{b_1}{\sqrt{2 p_1}})+\\ &&2\pi \left( T(\frac{n_2}{\sqrt{p_2}},\frac{\sqrt{p_2} o_1}{n_2}) + T(l_1,k_2)\right)\cdot erf(\frac{b_2}{\sqrt{2 p_2}})+\\ &&-2\pi\left( T(\frac{b_1}{\sqrt{p_1}},\frac{m_2}{\sqrt{\delta_2} b_1}) + T(\frac{b_2}{\sqrt{p_2}},\frac{m_1}{\sqrt{\delta_2} b_2})\right)+\\ &&4\pi \left( T^{(2)}(h,(a_j+\frac{b_j}{h})_{j=1}^2) -T^{(2)}(\frac{n_1}{\sqrt{p_1}},(\frac{b_1}{n_1},\frac{\sqrt{p_1} o_2}{n_1})) -T^{(2)}(\frac{n_2}{\sqrt{p_2}},(\frac{b_2}{n_2},\frac{\sqrt{p_2} o_1}{n_2})) -T^{(2)}(l_2,(\frac{\sqrt{\delta_2}b_1}{m_2},k_1)) -T^{(2)}(l_1,(\frac{\sqrt{\delta_2}b_2}{m_1},k_2)) \right) \end{eqnarray} As usual I include a piece of code that verifies this expression:

d = 2; M = 3;
mj =.; mj[i_] := If[i == 1, 2, 1];
Clear[phi]; phi[x_] := Exp[-1/2 x^2]/Sqrt[2 Pi];
Clear[T]; 
T[h_, a_, b_] := 
 1/(2 Pi) (ArcTan[a] - ArcTan[a + b/h] - ArcTan[a + h/b + a^2 h/b]) + 
  1/4 Erf[b/Sqrt[2 (1 + a^2)]] + OwenT[h, a + b/h] + 
  OwenT[b/Sqrt[1 + a^2], a + h/b + a^2 h/b];
gT[h_, a_] := 
  NIntegrate[
   phi[xi] Product[
     1/2 Erf[a[[i]]/Sqrt[2] xi], {i, 1, Length[a]}], {xi, h, 
    Infinity}];
For[count = 1, count <= 100, count++,
  h = RandomReal[{0, M}, WorkingPrecision -> 50];
  Clear[a]; 
  For[i = 1, i <= d, i++, 
   a[i] = RandomReal[{0, M}, WorkingPrecision -> 50]];
  Clear[b]; 
  For[i = 1, i <= d, i++, 
   b[i] = RandomReal[{0, M}, WorkingPrecision -> 50]];

  I1 = NIntegrate[
    phi[xi] Product[
      1/2 Erf[(a[i] xi + b[i])/Sqrt[2]], {i, 1, d}], {xi, h, 
     Infinity}];

  NIntegrate[
    phi[xi] Product[1/2 Erf[a[i]/Sqrt[2] xi], {i, 1, d}], {xi, 0, 
     Infinity}] - 
   NIntegrate[
    phi[xi] Product[
      1/2 Erf[(a[i] + 1/h b[i] )/Sqrt[2] xi], {i, 1, d}], {xi, 0, 
     Infinity}] + 
   NIntegrate[
    phi[xi] Product[
      1/2 Erf[(a[i] + 1/h b[i] )/Sqrt[2] xi], {i, 1, d}], {xi, h, 
     Infinity}] + 
   Sum[NIntegrate[
     phi[xi] T[(a[i] + h/b[i] (1 + a[i]^2)) xi, a[mj[i]]/Sqrt[
       1 + a[i]^2], (b[mj[i]] (1 + a[i]^2) - b[i] a[i] a[mj[i]])/(
        b[i] Sqrt[1 + a[i]^2]) xi], {xi, 0, b[i]/Sqrt[
      1 + a[i]^2]}], {i, 1, d}];
  dd = h^2 + (a[1] h + b[1])^2 + (a[2] h + b[2])^2;
  dd1 = h (1 + a[1]^2 + a[2]^2) + a[1] b[1] + a[2] b[2];
  dd2 = 1 + a[1]^2 + a[2]^2;


  {m1, m2} = {b[1] (1 + a[2]^2 ) - a[1] a[2] b[2], 
    b[2] (1 + a[1]^2 ) - a[1] a[2] b[1]};
  {n1, n2} = {h + h a[1]^2 + a[1] b[1], h + h a[2]^2 + a[2] b[2]};
  {o1, o2} = {h a[1] + b[1], h a[2] + b[2]};
  {p1, p2} = {1 + a[1]^2, 1 + a[2]^2};
  {k1, k2} = {(Sqrt[p1] (dd1))/m2, (Sqrt[p2] (dd1))/m1};
  {l1, l2} = {m1/Sqrt[(p2) (dd2)], m2/Sqrt[(p1) (dd2)]};


  I2 = 1/(
    4 \[Pi]) (ArcTan[(a[1] a[2])/Sqrt[dd2]] - 
      ArcTan[((o1) (o2))/(h Sqrt[dd])] +
      ArcTan[m2/(Sqrt[dd2] b[1])] + ArcTan[m1/(Sqrt[dd2] b[2])] +
      ArcTan[(b[1] (o2))/((n1) Sqrt[dd])] + 
      ArcTan[(b[2] (o1) )/((n2) Sqrt[dd])] +
      ArcTan[(b[1] (dd1))/Sqrt[(m2)^2 (dd)]] + 
      ArcTan[(b[2] (dd1))/Sqrt[(m1)^2 (dd)]] +
      (ArcTan[a[2]/Sqrt[p1]] - ArcTan[(Sqrt[p1] (o2))/n1] - 
         ArcTan[k1]) Erf[b[1]/(
        Sqrt[2] Sqrt[p1])] + (ArcTan[a[1]/Sqrt[p2]] - 
         ArcTan[(Sqrt[p2] (o1))/n2] - ArcTan[k2]) Erf[b[2]/(
        Sqrt[2] Sqrt[p2])] +
      2 \[Pi] (OwenT[n1/Sqrt[p1], (Sqrt[p1] (o2))/n1] + 
         OwenT[l2, k1]) Erf[b[1]/(Sqrt[2] Sqrt[p1])] +
      2 \[Pi] (OwenT[n2/Sqrt[p2], (Sqrt[p2] (o1))/n2] + 
         OwenT[l1, k2]) Erf[b[2]/(Sqrt[2] Sqrt[p2])] - 
      2 \[Pi] (OwenT[b[1]/Sqrt[p1], m2/(Sqrt[dd2] b[1])] + 
         OwenT[b[2]/Sqrt[p2], m1/(Sqrt[dd2] b[2])]) +
      4 \[Pi] (gT[h, {a[1] + b[1]/h, a[2] + b[2]/h}] +
         -gT[n1/Sqrt[p1], {b[1]/n1, (Sqrt[p1] (o2))/n1}] - 
         gT[n2/Sqrt[p2], {b[2]/n2, (Sqrt[p2] (o1))/n2}] - 
         gT[l2, {(Sqrt[dd2] b[1])/m2, k1}] - 
         gT[l1, {(Sqrt[dd2] b[2])/m1, k2}]));
  If[Abs[I2/I1 - 1] > 10^(-2), 
   Print["Results do not match..", {count, {a[1], a[2], b[1], b[2], 
      h}, {I1, I2}}]; Break[]];
  PrintTemporary[{count, I1, I2}];
  ];

Update: It might be interesting to know if it is possible to express the quantities $T^{(2)}(h,(a_1,a_2))$ in some alternative way. As a matter of fact by starting from the integral definition of this quantity , then differentiating with respect to $a_1$ and then integrating by parts and finally integrating with respect to $a_1$ from zero to $a_1$ we stumbled on a following formula: \begin{eqnarray} T^{(2)}(h,(a_1,a_2)) = \frac{2 \pi \text{erf}\left(\frac{\text{a2} h}{\sqrt{2}}\right) T(h,\text{a1})+\arctan\left(\frac{\text{a1} \text{a2}}{\sqrt{\text{a1}^2+\text{a2}^2+1}}\right) \text{erfc}\left(\frac{h \sqrt{\text{a1}^2+\text{a2}^2+1}}{\sqrt{2}}\right)}{4 \pi } + \frac{h \sqrt{1+a_2^2}}{\pi^{3/2} 2^{3/2}} \int\limits_0^{arccosh(\sqrt{\frac{1+a_1^2+a_2^2}{1+a_2^2}})} \sinh(\theta) \cdot \arctan\left(a_2 \frac{\sinh(\theta)}{\cosh(\theta)} \right) \cdot e^{-\frac{h^2}{2} (1+a_2^2) \cosh(\theta)^2} d\theta \end{eqnarray}

In particular for $h=0$ we have : \begin{equation} T^{(2)}(0,(a_1,a_2)) = \frac{1}{4 \pi} \arctan\left(\frac{\text{a1} \text{a2}}{\sqrt{\text{a1}^2+\text{a2}^2+1}}\right) \end{equation} as it should be (see An integral involving error functions and a Gaussian).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.