Multivariate gaussian integral over positive realsGaussian integrals over a half-spaceAn integral involving error functions and a GaussianGeneralized Owen's T functionAn integral involving a Gaussian, error functions and the Owen's T function.Derivation of a closed-form solution for the integral of a 3D Gaussian over the *positive* realsIntegration on a sphereTriple integral over an ellipsoidSurface area of $x^2+y^2+z^2=9$, where $1leq x^2+y^2leq4$ and $zgeq0$Evaluating this surface integralNumerical Multivariate Integration with respect to a Multivariate GaussianComponent square integration on the unit $n$-sphereCalculate the integral of $phi(x,y,z)$ over the surface of the area of the unit sphereIntegral of $prod_i=1^n (mathbfxcdot mathbfu_i) mathrmdmathbfx$ over the ($n-1$)-sphere.Evaluating a quartic integral that looks like expectation of multivariate normal random variableUnit vectors in vector integration of differential surfaces
How to split IPA spelling into syllables
Does capillary rise violate hydrostatic paradox?
Rendered textures different to 3D View
C++ lambda syntax
What properties make a magic weapon befit a Rogue more than a DEX-based Fighter?
Are hand made posters acceptable in Academia?
What is the purpose of using a decision tree?
Why does the frost depth increase when the surface temperature warms up?
Why is implicit conversion not ambiguous for non-primitive types?
Is divisi notation needed for brass or woodwind in an orchestra?
Air travel with refrigerated insulin
Asserting that Atheism and Theism are both faith based positions
Showing mass murder in a kid's book
How do you justify more code being written by following clean code practices?
How can I, as DM, avoid the Conga Line of Death occurring when implementing some form of flanking rule?
When is the exact date for EOL of Ubuntu 14.04 LTS?
Has the laser at Magurele, Romania reached a tenth of the Sun's power?
Would this string work as string?
Make a Bowl of Alphabet Soup
I keep switching characters, how do I stop?
Strange behavior in TikZ draw command
Taking the numerator and the denominator
Why is indicated airspeed rather than ground speed used during the takeoff roll?
What (if any) is the reason to buy in small local stores?
Multivariate gaussian integral over positive reals
Gaussian integrals over a half-spaceAn integral involving error functions and a GaussianGeneralized Owen's T functionAn integral involving a Gaussian, error functions and the Owen's T function.Derivation of a closed-form solution for the integral of a 3D Gaussian over the *positive* realsIntegration on a sphereTriple integral over an ellipsoidSurface area of $x^2+y^2+z^2=9$, where $1leq x^2+y^2leq4$ and $zgeq0$Evaluating this surface integralNumerical Multivariate Integration with respect to a Multivariate GaussianComponent square integration on the unit $n$-sphereCalculate the integral of $phi(x,y,z)$ over the surface of the area of the unit sphereIntegral of $prod_i=1^n (mathbfxcdot mathbfu_i) mathrmdmathbfx$ over the ($n-1$)-sphere.Evaluating a quartic integral that looks like expectation of multivariate normal random variableUnit vectors in vector integration of differential surfaces
$begingroup$
The multivariate gaussian integral over the whole $mathbfR^n$ has closed form solution
$$P = int_mathbfx in mathbfR^n exp left(-frac12 mathbfx^T mathbfA mathbfxright),dmathbfx = sqrtfrac(2pi)^ndet mathbfA$$
where $mathbfA$ is a symmetric positive-definite covariance matrix.
However, I need to solve the integral for positive reals $mathbfx in mathbfR^n :, mathbfx_i geq 0 forall i$ only and in at least 6 dimensions:
$$P = int_mathbfx in mathbfR^n :, mathbfx_i geq 0 forall i exp left(-frac12 mathbfx^T mathbfA mathbfxright),dmathbfx$$
For diagonal $mathbfA$ with zero covariance, a solution has been published.
For non-diagonal covariance, my approach so far is to apply affine coordinate transforms to rotate and rescale the gaussian ellipsoid into the unit sphere (see here).
In two dimensions, the solution to the integral then reduces to comparing the area enclosed by the transformed positive coordinate axes (blue) to the area of the unit circle:
In three dimensions, the solution is given by the ration of the surface area of an enclosed spherical polygon to the surface area of the unit sphere.
In four dimensions, this approach becomes quite complicated, and I don't know how to use the usual spherical excess formulas for higher dimensions.
Any ideas or alternative approaches? Is there a multivariate error function? Any treatment on the multivariate half normal distribution?
integration
edited Apr 13 '17 at 12:58
Community♦
1
asked Jul 17 '14 at 1:05
le_mle_m
15816
$endgroup$
add a comment |
$begingroup$
The multivariate gaussian integral over the whole $mathbfR^n$ has closed form solution
$$P = int_mathbfx in mathbfR^n exp left(-frac12 mathbfx^T mathbfA mathbfxright),dmathbfx = sqrtfrac(2pi)^ndet mathbfA$$
where $mathbfA$ is a symmetric positive-definite covariance matrix.
However, I need to solve the integral for positive reals $mathbfx in mathbfR^n :, mathbfx_i geq 0 forall i$ only and in at least 6 dimensions:
$$P = int_mathbfx in mathbfR^n :, mathbfx_i geq 0 forall i exp left(-frac12 mathbfx^T mathbfA mathbfxright),dmathbfx$$
For diagonal $mathbfA$ with zero covariance, a solution has been published.
For non-diagonal covariance, my approach so far is to apply affine coordinate transforms to rotate and rescale the gaussian ellipsoid into the unit sphere (see here).
In two dimensions, the solution to the integral then reduces to comparing the area enclosed by the transformed positive coordinate axes (blue) to the area of the unit circle:
In three dimensions, the solution is given by the ration of the surface area of an enclosed spherical polygon to the surface area of the unit sphere.
In four dimensions, this approach becomes quite complicated, and I don't know how to use the usual spherical excess formulas for higher dimensions.
Any ideas or alternative approaches? Is there a multivariate error function? Any treatment on the multivariate half normal distribution?
integration
edited Apr 13 '17 at 12:58
Community♦
1
asked Jul 17 '14 at 1:05
le_mle_m
15816
$endgroup$
$begingroup$
Your notation, $mathbf xgeq0$ doesn't make sense when $mathbfxinmathbbR^n$, $ngeq 2$.
$endgroup$
– nullgeppetto
Jul 17 '14 at 14:21
$begingroup$
I think that you should use the solution provided in the link you give about my older question.
$endgroup$
– nullgeppetto
Jul 17 '14 at 14:23
$begingroup$
With $x geq 0$ I meant $x_i geq 0$ $forall i$. Is this bad notation? Your question is different from mine in that your integral can be reduced to n-1 gaussian integrals and one "half gaussian" one, which can be solved by the error function. I am not aware of an n-dimensional error function though...
$endgroup$
– le_m
Jul 17 '14 at 14:32
$begingroup$
I encountered similar problem recently. For the case of all covariances being the same, the problem is exactly solvable. If this is the one you are interested in, let me know.
$endgroup$
– Sungmin
Nov 6 '15 at 14:55
$begingroup$
@Sungmin: Would you mind showing me the solution in this case ? Honestly I truly doubt that any results exist when $n > 4$. For $n=4$ di-logarithms enter the result and in case of higher dimensions higher order poly-loarithms will be needed.
$endgroup$
– Przemo
Mar 12 at 18:35
add a comment |
$begingroup$
The multivariate gaussian integral over the whole $mathbfR^n$ has closed form solution
$$P = int_mathbfx in mathbfR^n exp left(-frac12 mathbfx^T mathbfA mathbfxright),dmathbfx = sqrtfrac(2pi)^ndet mathbfA$$
where $mathbfA$ is a symmetric positive-definite covariance matrix.
However, I need to solve the integral for positive reals $mathbfx in mathbfR^n :, mathbfx_i geq 0 forall i$ only and in at least 6 dimensions:
$$P = int_mathbfx in mathbfR^n :, mathbfx_i geq 0 forall i exp left(-frac12 mathbfx^T mathbfA mathbfxright),dmathbfx$$
For diagonal $mathbfA$ with zero covariance, a solution has been published.
For non-diagonal covariance, my approach so far is to apply affine coordinate transforms to rotate and rescale the gaussian ellipsoid into the unit sphere (see here).
In two dimensions, the solution to the integral then reduces to comparing the area enclosed by the transformed positive coordinate axes (blue) to the area of the unit circle:
In three dimensions, the solution is given by the ration of the surface area of an enclosed spherical polygon to the surface area of the unit sphere.
In four dimensions, this approach becomes quite complicated, and I don't know how to use the usual spherical excess formulas for higher dimensions.
Any ideas or alternative approaches? Is there a multivariate error function? Any treatment on the multivariate half normal distribution?
integration
edited Apr 13 '17 at 12:58
Community♦
1
asked Jul 17 '14 at 1:05
le_mle_m
15816
$endgroup$
The multivariate gaussian integral over the whole $mathbfR^n$ has closed form solution
$$P = int_mathbfx in mathbfR^n exp left(-frac12 mathbfx^T mathbfA mathbfxright),dmathbfx = sqrtfrac(2pi)^ndet mathbfA$$
where $mathbfA$ is a symmetric positive-definite covariance matrix.
However, I need to solve the integral for positive reals $mathbfx in mathbfR^n :, mathbfx_i geq 0 forall i$ only and in at least 6 dimensions:
$$P = int_mathbfx in mathbfR^n :, mathbfx_i geq 0 forall i exp left(-frac12 mathbfx^T mathbfA mathbfxright),dmathbfx$$
For diagonal $mathbfA$ with zero covariance, a solution has been published.
For non-diagonal covariance, my approach so far is to apply affine coordinate transforms to rotate and rescale the gaussian ellipsoid into the unit sphere (see here).
In two dimensions, the solution to the integral then reduces to comparing the area enclosed by the transformed positive coordinate axes (blue) to the area of the unit circle:
In three dimensions, the solution is given by the ration of the surface area of an enclosed spherical polygon to the surface area of the unit sphere.
In four dimensions, this approach becomes quite complicated, and I don't know how to use the usual spherical excess formulas for higher dimensions.
Any ideas or alternative approaches? Is there a multivariate error function? Any treatment on the multivariate half normal distribution?
integration
integration
edited Apr 13 '17 at 12:58
Community♦
1
asked Jul 17 '14 at 1:05
le_mle_m
15816
edited Apr 13 '17 at 12:58
Community♦
1
asked Jul 17 '14 at 1:05
le_mle_m
15816
edited Apr 13 '17 at 12:58
Community♦
1
edited Apr 13 '17 at 12:58
Community♦
1
edited Apr 13 '17 at 12:58
Community♦
1
1
asked Jul 17 '14 at 1:05
le_mle_m
15816
asked Jul 17 '14 at 1:05
le_mle_m
15816
asked Jul 17 '14 at 1:05
le_mle_m
15816
15816
$begingroup$
Your notation, $mathbf xgeq0$ doesn't make sense when $mathbfxinmathbbR^n$, $ngeq 2$.
$endgroup$
– nullgeppetto
Jul 17 '14 at 14:21
$begingroup$
I think that you should use the solution provided in the link you give about my older question.
$endgroup$
– nullgeppetto
Jul 17 '14 at 14:23
$begingroup$
With $x geq 0$ I meant $x_i geq 0$ $forall i$. Is this bad notation? Your question is different from mine in that your integral can be reduced to n-1 gaussian integrals and one "half gaussian" one, which can be solved by the error function. I am not aware of an n-dimensional error function though...
$endgroup$
– le_m
Jul 17 '14 at 14:32
$begingroup$
I encountered similar problem recently. For the case of all covariances being the same, the problem is exactly solvable. If this is the one you are interested in, let me know.
$endgroup$
– Sungmin
Nov 6 '15 at 14:55
$begingroup$
@Sungmin: Would you mind showing me the solution in this case ? Honestly I truly doubt that any results exist when $n > 4$. For $n=4$ di-logarithms enter the result and in case of higher dimensions higher order poly-loarithms will be needed.
$endgroup$
– Przemo
Mar 12 at 18:35
add a comment |
$begingroup$
Your notation, $mathbf xgeq0$ doesn't make sense when $mathbfxinmathbbR^n$, $ngeq 2$.
$endgroup$
– nullgeppetto
Jul 17 '14 at 14:21
$begingroup$
I think that you should use the solution provided in the link you give about my older question.
$endgroup$
– nullgeppetto
Jul 17 '14 at 14:23
$begingroup$
With $x geq 0$ I meant $x_i geq 0$ $forall i$. Is this bad notation? Your question is different from mine in that your integral can be reduced to n-1 gaussian integrals and one "half gaussian" one, which can be solved by the error function. I am not aware of an n-dimensional error function though...
$endgroup$
– le_m
Jul 17 '14 at 14:32
$begingroup$
I encountered similar problem recently. For the case of all covariances being the same, the problem is exactly solvable. If this is the one you are interested in, let me know.
$endgroup$
– Sungmin
Nov 6 '15 at 14:55
$begingroup$
@Sungmin: Would you mind showing me the solution in this case ? Honestly I truly doubt that any results exist when $n > 4$. For $n=4$ di-logarithms enter the result and in case of higher dimensions higher order poly-loarithms will be needed.
$endgroup$
– Przemo
Mar 12 at 18:35
$begingroup$
Your notation, $mathbf xgeq0$ doesn't make sense when $mathbfxinmathbbR^n$, $ngeq 2$.
$endgroup$
– nullgeppetto
Jul 17 '14 at 14:21
$begingroup$
Your notation, $mathbf xgeq0$ doesn't make sense when $mathbfxinmathbbR^n$, $ngeq 2$.
$endgroup$
– nullgeppetto
Jul 17 '14 at 14:21
$begingroup$
I think that you should use the solution provided in the link you give about my older question.
$endgroup$
– nullgeppetto
Jul 17 '14 at 14:23
$begingroup$
I think that you should use the solution provided in the link you give about my older question.
$endgroup$
– nullgeppetto
Jul 17 '14 at 14:23
$begingroup$
With $x geq 0$ I meant $x_i geq 0$ $forall i$. Is this bad notation? Your question is different from mine in that your integral can be reduced to n-1 gaussian integrals and one "half gaussian" one, which can be solved by the error function. I am not aware of an n-dimensional error function though...
$endgroup$
– le_m
Jul 17 '14 at 14:32
$begingroup$
With $x geq 0$ I meant $x_i geq 0$ $forall i$. Is this bad notation? Your question is different from mine in that your integral can be reduced to n-1 gaussian integrals and one "half gaussian" one, which can be solved by the error function. I am not aware of an n-dimensional error function though...
$endgroup$
– le_m
Jul 17 '14 at 14:32
$begingroup$
I encountered similar problem recently. For the case of all covariances being the same, the problem is exactly solvable. If this is the one you are interested in, let me know.
$endgroup$
– Sungmin
Nov 6 '15 at 14:55
$begingroup$
I encountered similar problem recently. For the case of all covariances being the same, the problem is exactly solvable. If this is the one you are interested in, let me know.
$endgroup$
– Sungmin
Nov 6 '15 at 14:55
$begingroup$
@Sungmin: Would you mind showing me the solution in this case ? Honestly I truly doubt that any results exist when $n > 4$. For $n=4$ di-logarithms enter the result and in case of higher dimensions higher order poly-loarithms will be needed.
$endgroup$
– Przemo
Mar 12 at 18:35
$begingroup$
@Sungmin: Would you mind showing me the solution in this case ? Honestly I truly doubt that any results exist when $n > 4$. For $n=4$ di-logarithms enter the result and in case of higher dimensions higher order poly-loarithms will be needed.
$endgroup$
– Przemo
Mar 12 at 18:35
add a comment |
5 Answers
5
active
oldest
votes
$begingroup$
The integral over (coordinate-wise) positive values appears in the treatment of dichotomized Gaussian distributions, so you might find the answer to your problem there. Relevant references would be:
- DR Cox, N Wermuth, Biometrika, 2002
- JH Macke, P Berens et al., Neural Computation, 2009
answered Jul 17 '14 at 14:51
flipsflips
515
$endgroup$
$begingroup$
Thanks! I looked into your referenced papers and it seems they approximate the integral in question under various assumptions, e.g. all covariances being identical.
$endgroup$
– le_m
Aug 4 '14 at 23:28
$begingroup$
As far as I can see, these only cite a solution for the case n=2, tracing back to Sheppard 1898
$endgroup$
– guillefix
Feb 13 at 15:09
add a comment |
$begingroup$
Let us compute the result in case $n=2$. Here the matrix reads $A=left(beginarrayrra & c\c& bendarrayright)$ .Therefore we have:
begineqnarray
P&=& intlimits_mathbb R_+^2 expleft-frac12left[sqrta(s_1+fracca s_2)right]^2 -frac12 fracb a-c^2a s_2^2right ds_1 ds_2\
&=&frac1sqrta sqrtfracpi2 intlimits_0^infty erfcleft(fraccsqrta fracs_2sqrt2 right)expleft-frac12(fracb a-c^2a)s_2^2 rightds_2\
&=&sqrtfracpi2 frac1sqrtb a-c^2 intlimits_0^infty erfc(fraccsqrtb a-c^2 fracs_2sqrt2) e^-frac12 s_2^2 ds_2\
&=& sqrtfracpi2 frac1sqrtb a-c^2 left(
sqrtfracpi2- sqrtfrac2pi arctan(fraccsqrtb a-c^2)right)\
&=& frac1sqrtb a-c^2 arctan(fracsqrtb a-c^2c)
endeqnarray
In the top line we completed the first integration variable to a square and in the second line we integrated over that variable. In the third line we changed variables accordingly . In the fourth line we integrated over the second variable by writing $erfc() = 1- erf()$ and then expanding the error function in a Taylor series and integrating term by term and finally in the last line we simplified the result.
Now, by doing similar calculations we obtained the following result in case $n=3$. Here $A=left(beginarrayrrra & a_12 & a_13\a_12& b&a_23\a_13&a_23&cendarrayright)$.
Firstly we have:
begineqnarray
&&vecs^(T).(A.vecs) = \
&&left(sqrta ( s_1 + fraca_1,2 s_2 + a_1,3 s_3a )right)^2 +
left( b- fraca_1,2^2aright) s_2^2 + left(c-fraca_1,3^2aright) s_3^2 + 2 left(a_2,3-fraca_1,2 a_1,3 aright) s_2 s_3
endeqnarray
Therefore integrating over $s_1$ gives:
begineqnarray
&&P=sqrtfracpi 2 frac1sqrta cdot \
&&intlimits_bf R^2 texterfcleft(fraca_1,2 s_2+a_1,3 s_3sqrt2 sqrtaright) cdot \
&&exp left[
-frac12 left(s_2^2 left(b-fraca_1,2^2aright)+2
s_2 s_3 left(a_2,3-fraca_1,2 a_1,3aright)+s_3^2 left(c-fraca_1,3^2aright)right)
right]
ds_2 ds_3=\
&&
fracsqrtpi a_1,2 intlimits_0^infty
texterfc(u) cdot expleft[-frac12u^2 (frac2 a ba_1,2^2 - 2)right]\
&&
intlimits_0^fracsqrt2 aa_1,3 u
exp left[-frac12 left(s_3 ufrac2 sqrt2 sqrta a_1,2 left(a_2,3-fracb
a_1,3a_1,2right)+
s_3^2fraca_1,3 a_1,2 left(fraca_1,3 ba_1,2+fraca_1,2 ca_1,3-2 a_2,3right)right)right] ds_3 du
endeqnarray
Now it is clear that we can do the integral over $s_3$ in the sense that we can express it through a difference of error functions.Denote $delta:=-2 a_1,2 a_1,3 a_2,3 +a_1,3^2 b +a_1,2^2 c$. Then we have
begineqnarray
&&P=fracpisqrt2sqrtdelta cdotintlimits_0^infty erfc(u) left( erfleft[fracsqrta(-a_1,3 a_2,3+a_1,2 c)a_1,3 sqrtdelta u right] - erfleft[ fracsqrta(a_1,2 a_2,3-a_1,3 b)a_1,2 sqrtdelta u right]right) e^-fracdet(A) delta u^2 du=\
&&fracpisqrt2 det(A)cdot \
&&
intlimits_0^infty erfcleft(u sqrtfracdeltadet(A)right)e^-u^2cdot \
&&left(-erfc(sqrta frac(-a_13a_23+a_12 c)a_13 sqrtdet(A) u)+erfc(sqrta frac(a_12a_23-a_13 b)a_12 sqrtdet(A) u)right) du \
&&=sqrtfracpi2 det(A)\
left[right.\
&&-arctanleft(fraca_13 sqrtdet(A)sqrta(-a_13a_23+a_12 c)right)+
arctanleft(fracsqrtc sqrtdet(A)-a_13 a_23 + a_12 cright)
\
&&+arctanleft(fraca_12 sqrtdet(A)sqrta (a_12 a_23 - a_13 b)right)-arctanleft(fracsqrtb sqrtdet(A)a_12 a_23 - a_13 bright)
left. right]\
&&=sqrtfracpi2 det(A)\
&&left[right.\
&&left.
arctanleft(frac(a_1,3-sqrta_1,1a_3,3)(a_1,3a_2,3-a_1,2a_3,3)sqrta_1,1 (a_1,3a_2,3-a_1,2a_3,3)^2+a_1,3 sqrta_3,3 det(A) sqrtdet(A)right)+right.\
&&left.
arctanleft(frac(a_1,2-sqrta_1,1a_2,2)(a_1,2a_2,3-a_1,3a_2,2)sqrta_1,1 (a_1,2a_2,3-a_1,3a_2,2)^2+a_1,2 sqrta_2,2 det(A) sqrtdet(A)right)
right]
endeqnarray
where in the last line we used An integral involving error functions and a Gaussian .
I also include a Mathematica code snippet that verifies all the steps involved:
(*3d*)
A =.; B =.; CC =.; A12 =.; A23 =.; A13 =.;
For[DDet = 0, True, ,
A, B, CC, A12, A23, A13 =
RandomReal[0, 1, 6, WorkingPrecision -> 50];
DDet = Det[A, A12, A13, A12, B, A23, A13, A23, CC];
If[DDet > 0, Break[]];
];
a = Sqrt[(-2 A12 A13 A23 + A13^2 B + A12^2 CC)/DDet];
b1, b2 = ( Sqrt[A] (-A13 A23 + A12 CC))/ Sqrt[DDet], (
Sqrt[A] (A12 A23 - A13 B))/ Sqrt[DDet];
AA1, AA2 = 2 Sqrt[2] Sqrt[
A] (( A23 A12 - A13 B)/A12^2), (-2 A12 A13 A23 + A13^2 B +
A12^2 CC)/A12^2;
DDet, a, b1, b2;
NIntegrate[
Exp[-1/2 (A s1^2 + B s2^2 + CC s3^2 + 2 A12 s1 s2 + 2 A23 s2 s3 +
2 A13 s1 s3)], s1, 0, Infinity, s2, 0, Infinity, s3, 0,
Infinity]
NIntegrate[
Exp[-1/2 ((Sqrt[A] (s1 + (A12 s2 + A13 s3)/A))^2 + (B -
A12^2/A) s2^2 + (CC - A13^2/A) s3^2 +
2 (A23 - A12 A13/A) s2 s3)], s1, 0, Infinity, s2, 0,
Infinity, s3, 0, Infinity]
NIntegrate[
1/Sqrt[A] Sqrt[
Pi/2] Erfc[(A12 s2 + A13 s3)/
Sqrt[2 A]] Exp[-1/
2 ((B - A12^2/A) s2^2 + (CC - A13^2/A) s3^2 +
2 (A23 - A12 A13/A) s2 s3)], s2, 0, Infinity, s3, 0,
Infinity]
Sqrt[Pi]/A12 NIntegrate[
Erfc[u] Exp[-1/
2 ( A13/A12 (-2 A23 + (A13 B)/A12 + CC A12/A13) s3^2 + (
2 Sqrt[2] Sqrt[A] )/
A12 ( A23 - ( A13 B)/A12) s3 u + (-2 + (2 A B)/
A12^2) u^2)], u, 0, Infinity, s3, 0, Sqrt[2 A]/A13 u]
Sqrt[Pi]/A12 NIntegrate[
Erfc[u] Exp[-1/2 (Sqrt[AA2] s3 + u/2 AA1/Sqrt[AA2])^2] Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity, s3, 0, Sqrt[2 A]/A13 u]
Sqrt[Pi]/(A12 Sqrt[AA2])
NIntegrate[
Erfc[u] Exp[-1/2 (s3)^2] Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity, s3,
u/2 AA1/Sqrt[AA2], ((A13 AA1 + 2 AA2 Sqrt[2] Sqrt[A]) u)/(
2 A13 Sqrt[AA2])]
Sqrt[Pi]/(A12 Sqrt[AA2]) Sqrt[[Pi]/2]
NIntegrate[
Erfc[u] (
Erf[(A13 AA1 + 2 AA2 Sqrt[2] Sqrt[A])/(2 A13 Sqrt[2] Sqrt[AA2])
u] - Erf[AA1/(2 Sqrt[2] Sqrt[AA2]) u]) Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity]
Pi/Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC] Sqrt[1/2]
NIntegrate[
Erfc[u] (
Erf[( Sqrt[A] (-A13 A23 + A12 CC) u)/(
A13 Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC])] -
Erf[(Sqrt[A] (A12 A23 - A13 B) u)/(
A12 Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC])]) Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity]
Pi/ Sqrt[-2 A12 A13 A23 + A13^2 B +
A12^2 CC] Sqrt[1/2] a NIntegrate[
Erfc[a u] (
Erf[( Sqrt[A] (-A13 A23 + A12 CC) u)/(A13 Sqrt[DDet])] -
Erf[(Sqrt[A] (A12 A23 - A13 B) u)/(A12 Sqrt[DDet])]) Exp[-
u^2], u, 0, Infinity]
Pi/Sqrt[2 DDet] NIntegrate[(Erfc[u a]) Exp[-u^2] (Erf[b1/A13 u] -
Erf[b2/A12 u]), u, 0, Infinity]
Sqrt[Pi]/Sqrt[
2 DDet] (ArcTan[ Sqrt[A]/A13 (-A13 A23 + A12 CC)/ Sqrt[DDet]] -
ArcTan[1/ Sqrt[CC] (-A13 A23 + A12 CC)/ Sqrt[DDet]] -
ArcTan[ Sqrt[A]/A12 (A12 A23 - A13 B)/ Sqrt[DDet]] +
ArcTan[ 1/Sqrt[B] (A12 A23 - A13 B)/ Sqrt[DDet]])
-(Sqrt[Pi]/
Sqrt[2 DDet]) (ArcTan[(A13 Sqrt[DDet])/(
Sqrt[A] (-A13 A23 + A12 CC))] -
ArcTan[(Sqrt[CC] Sqrt[DDet])/(-A13 A23 + A12 CC)] -
ArcTan[(A12 Sqrt[DDet])/(Sqrt[A] (A12 A23 - A13 B))] +
ArcTan[(Sqrt[B] Sqrt[DDet])/(A12 A23 - A13 B)])
Sqrt[Pi]/Sqrt[
2 DDet] (ArcTan[((A13 - Sqrt[A] Sqrt[CC]) (A13 A23 - A12 CC) Sqrt[
DDet])/(Sqrt[A] (A13 A23 - A12 CC)^2 + A13 Sqrt[CC] DDet)] +
ArcTan[((A12 - Sqrt[A] Sqrt[B]) (A12 A23 - A13 B) Sqrt[DDet])/(
Sqrt[A] (A12 A23 - A13 B)^2 + A12 Sqrt[B] DDet)])
Update: Now let us take a look at the $n=4$ case. In here:
beginequation
bf A=left(
beginarrayrrrr
a & a_1,2 & a_1,3 & a_1,4 \
a_1,2 & b & a_2,3 & a_2,4 \
a_1,3 & a_2,3 & c & a_3,4 \
a_1,4 & a_2,4 & a_3,4 & d
endarray
right)
endequation
then by doing basically the same computations as a above we managed to reduce the integral in question to a following two dimensional integral. We have:
begineqnarray
&&P= \
&&!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!fracpisqrt2 delta intlimits_0^infty intlimits_0^fracsqrt2 aa_1,2 u
erfc[u] cdot expleft[fracmathfrak A_0,0 u^2 + mathfrak A_1,0 u s_2 +mathfrak A_1,1 s_2^22 deltaright] cdot left( erf[fracmathfrak B_1 u + mathfrak B_2 s_2a_1,3 sqrt2 delta] + erf[fracmathfrak C_1 u + mathfrak C_2 s_2a_1,4 sqrt2 delta]right) d s_2 du=\
&&frac2 imath pi^3/2sqrtmathfrak A_1,1
intlimits_0^infty erfc[u] expfrac4 mathfrak A_0,0 mathfrak A_1,1 - mathfrak A_1,0^28delta mathfrak A_1,1 u^2cdot\
&&
left[right.\
&&left.
left.Tleft(frac(mathfrak A_1,0+xi) u2imath sqrtmathfrak A_1,1 delta,
fracimath mathfrak B_2a_1,3 sqrtmathfrak A_1,1,fracu(2mathfrak A_1,1 mathfrak B_1 - mathfrak A_1,0 mathfrak B_2)2sqrtdelta a_1,3 mathfrak A_1,1right)right|_frac2mathfrak A_1,1 sqrt2 aa_1,2^0 +right.\
&&left.
left.Tleft(frac(mathfrak A_1,0+xi) u2imath sqrtmathfrak A_1,1 delta,
fracimath mathfrak C_2a_1,3 sqrtmathfrak A_1,1,fracu(2mathfrak A_1,1 mathfrak C_1 - mathfrak A_1,0 mathfrak C_2)2sqrtdelta a_1,3 mathfrak A_1,1right)right|_frac2mathfrak A_1,1 sqrt2 aa_1,2^0 +right.\
&&left.
right] du quad (i)
endeqnarray
where $T(cdot,cdot,cdot)$ is the generalized Owen's T function Generalized Owen's T function and
begineqnarray
delta&:=&a_1,3(a_1,3 d-a_1,4 a_3,4) + a_1,4(a_1,4 c- a_1,3 a_3,4)\
mathfrak A_0,0&:=&2 a left(a_3,4^2-c dright)+2 a_1,4 (a_1,4 c-a_1,3 a_3,4)+2 a_1,3 (a_1,3 d-a_1,4 a_3,4)\
mathfrak A_1,0&:=&2 sqrt2 sqrta left(a_1,2 left(c d-a_3,4^2right)+a_1,3
(a_2,4 a_3,4-a_2,3 d)+a_1,4 (a_2,3 a_3,4-a_2,4 c)right)\
mathfrak A_1,1&:=&a_1,2^2 left(a_3,4^2-c dright)+2 a_1,2 a_1,3 (a_2,3 d-a_2,4 a_3,4)+2 a_1,2 a_1,4
(a_2,4 c-a_2,3 a_3,4)+a_1,3^2 left(a_2,4^2-b dright)+2 a_1,3 a_1,4 (a_3,4 b-a_2,3 a_2,4)+a_1,4^2 left(a_2,3^2-b cright)\
hline\
mathfrak B_1&:=&sqrt2 sqrta (a_1,4
c-a_1,3 a_3,4)\
mathfrak B_2&:=&a_1,2 (a_1,3 a_3,4-a_1,4 c)+a_1,3 (a_1,4 a_2,3-a_1,3 a_2,4)\
mathfrak C_1&:=&sqrt2 sqrta (a_1,3 d-a_1,4 a_3,4)\
mathfrak C_2&:=&a_1,2 (a_1,4
a_3,4-a_1,3 d)+a_1,4 (a_1,3 a_2,4-a_1,4 a_2,3)
endeqnarray
nu = 4; Clear[T]; Clear[a]; x =.;
(*a0.dat, a1.dat or a2.dat*)
mat = << "a0.dat";
a, b, c, d, a12, a13, a14, a23, a24, a34 = mat[[1, 1]],
mat[[2, 2]], mat[[3, 3]], mat[[4, 4]], mat[[1, 2]], mat[[1, 3]],
mat[[1, 4]], mat[[2, 3]], mat[[2, 4]], mat[[3, 4]];
dd, A00, A10,
A11 = -2 a13 a14 a34 + a14^2 c + a13^2 d, -4 a13 a14 a34 +
2 a a34^2 + 2 a14^2 c + 2 a13^2 d - 2 a c d,
2 Sqrt[2] Sqrt[a] a14 a23 a34 + 2 Sqrt[2] Sqrt[a] a13 a24 a34 -
2 Sqrt[2] Sqrt[a] a12 a34^2 - 2 Sqrt[2] Sqrt[a] a14 a24 c -
2 Sqrt[2] Sqrt[a] a13 a23 d + 2 Sqrt[2] Sqrt[a] a12 c d,
a14^2 a23^2 - 2 a13 a14 a23 a24 + a13^2 a24^2 -
2 a12 a14 a23 a34 - 2 a12 a13 a24 a34 + a12^2 a34^2 +
2 a13 a14 a34 b + 2 a12 a14 a24 c - a14^2 b c + 2 a12 a13 a23 d -
a13^2 b d - a12^2 c d;
B1, B2, C1,
C2 = Sqrt[2] Sqrt[
a] (-a13 a34 + a14 c), (a13 a14 a23 - a13^2 a24 + a12 a13 a34 -
a12 a14 c),
Sqrt[2] Sqrt[
a] (-a14 a34 + a13 d), (-a14^2 a23 + a13 a14 a24 + a12 a14 a34 -
a12 a13 d);
NIntegrate[
Exp[-1/2 Sum[mat[[i, j]] s[i] s[j], i, 1, nu, j, 1, nu]],
Evaluate[Sequence @@ Table[s[eta], 0, Infinity, eta, 1, nu]]]
Sqrt[[Pi]/(2 a)]
NIntegrate[
Erfc[(a12 s[2] + a13 s[3] + a14 s[4])/Sqrt[
2 a]] Exp[-1/
2 ((-(a12^2/a) + b) s[2]^2 + (-(a13^2/a) + c) s[
3]^2 + (-(a14^2/a) + d) s[4]^2 +
2 (-(( a13 a14)/a) + a34) s[3] s[4] +
2 (-(( a12 a13)/a) + a23) s[2] s[3] +
2 (-(( a12 a14)/a) + a24) s[2] s[4])],
Evaluate[Sequence @@ Table[s[eta], 0, Infinity, eta, 2, nu]]]
Sqrt[[Pi]]
1/a14 NIntegrate[
Erfc[u] Exp[(
2 a14 a24 s[2] (-Sqrt[2] Sqrt[a] u + a12 s[2]) -
d (2 a u^2 - 2 Sqrt[2] Sqrt[a] a12 u s[2] + a12^2 s[2]^2) +
a14^2 (2 u^2 - b s[2]^2))/(
2 a14^2) + ((Sqrt[2] Sqrt[
a] (-a14 a34 + a13 d) u + (-a14^2 a23 + a13 a14 a24 +
a12 a14 a34 - a12 a13 d) s[2]) s[3])/
a14^2 - ((-2 a13 a14 a34 + a14^2 c + a13^2 d) s[3]^2)/(
2 a14^2)], u, 0, Infinity, s[2], 0,
Sqrt[2] Sqrt[a]/a12 u, s[3], 0, (Sqrt[2 a] u - a12 s[2])/a13]
Pi/Sqrt[2 dd]
NIntegrate[
Erfc[u] Exp[(A00 u^2 + A10 u s[2] + A11 s[2]^2)/(
2 (dd))] (Erf[(B1 u + B2 s[2])/( a13 Sqrt[2 dd])] +
Erf[(C1 u + C2 s[2])/( a14^1 Sqrt[2 dd])]), u, 0,
Infinity, s[2], 0, Sqrt[2] Sqrt[a]/a12 u]
Now, I will provide the result. Note that the only assumptions on the underlying matrix $bf A$ are that it is symmetric and that its elements are non-negative. Firstly let us define:
begineqnarray
&&mathfrak J^(1,1)(a,b,c)= frac1pi^2cdot left(right.\
&&left.
-frac18
sumlimits_i=1^4 sumlimits_j=1^4
(-1)^j-1+lfloor fraci-12 rfloor
%
mathfrak F^(1,fracsqrt1+2 a^2+b^2 - sqrt2 asqrt1+b^2)_fraci sqrtb^2 c^2+b^2+1 (-1)^leftlfloor
fracj-12rightrfloor +i b c (-1)^jsqrtb^2+1,-fracb
(-1)^i+i (-1)^leftlceil fraci-12rightrceil
sqrtb^2+1
%
right. \
&&left.
right)quad (ii)
endeqnarray
where $mathfrak F^(A,B)_a,b$ is related to di-logarithms and is defined in An integral involving a Gaussian, error functions and the Owen's T function. .
Then we define another function as follows:
beginequation
bar mathfrak J^(1,1)(a,b,c):= fracpi2 arctanleft[ fracsqrt2 a csqrt2 a+b^2(1+c^2)right] - fracpi2 arctanleft[ cright] - 2 pi^2 mathfrak J^(1,1)(frac1sqrt2 a,fracbsqrt2 a,c)
endequation
and then the following quantities that depend on the underlying matrix. We have:
begineqnarray
delta&:=& a_3,3 a_4,1^2 - 2 a_3,1 a_3,4 a_4,1 + a_4,4 a_3,1^2\
W&:=&left(a_3,3 a_4,4-a_3,4^2right) a_1,2^2+2 a_1,4 (a_2,3
a_3,4-a_2,4 a_3,3) a_1,2+2 a_1,3 (a_2,4 a_3,4-a_2,3 a_4,4)
a_1,2+a_1,4^2 left(a_2,2 a_3,3-a_2,3^2right)+2 a_1,3 a_1,4
(a_2,3 a_2,4-a_2,2 a_3,4)+a_1,3^2 left(a_2,2
a_4,4-a_2,4^2right)\
W_1&:=&2 sqrta_1,1 left(a_1,4 (a_2,4 a_3,3-a_2,3 a_3,4)+a_1,3 (a_2,3
a_4,4-a_2,4 a_3,4)+a_1,2 left(a_3,4^2-a_3,3 a_4,4right)right)\
%
v_1&:=&frac1a_4,1 sqrtdelta left( sqrta_1,1(a_3,4 a_4,1 - a_3,1 a_4,4),-a_2,4 a_3,1 a_4,1 + a_2,3 a_4,1^2+a_2,1(-a_3,4 a_4,1+a_3,1a_4,4)right)\
v_2&:=&-frac1a_3,1 sqrtdelta left(sqrta_1,1(a_3,4 a_3,1 - a_4,1 a_3,3),-a_3,1 a_3,2 a_4,1 +a_2,4 a_3,1^2 + a_2,1(-a_3,4 a_3,1+a_4,1a_3,3) right)\
%
left( A, B right)&:=& frac1delta left( W,W_1 right)\
left( bf a_1,bf a_2 right)&:=& frac1sqrtA left(v_1(2),v_2(2) right)\
bf b_1&:=& sqrt2 v_1(1) - fracBsqrt2 A v_1(2)\
bf b_2&:=& sqrt2 v_2(1) - fracBsqrt2 A v_2(2)\
x&:=& fracsqrta_1,1a_2,1
endeqnarray
Then the result reads:
begineqnarray
&&P=frac1det(bf A) left(right.\
%
&& bar mathfrak J^(1,1)left( fracdet(bf A)W,fracBsqrt2 A,bf a_2+fracsqrt2 A bf b_2Bright) -
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracB+2 A xsqrt2 A,bf a_2+fracsqrt2 A bf b_2B+2 A xright)+\
&&!!!!!!!!!!!!!!!!!!!!!!!!!!! bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_2sqrt1+bf a_2^2,bf a_2+fracB(1+bf a_2^2)sqrt2 Abf b_2right) -
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_2sqrt1+bf a_2^2,bf a_2+frac(B+2 A x)(1+bf a_2^2)sqrt2 Abf b_2right)+\
%
&& -bar mathfrak J^(1,1)left( fracdet(bf A)W,fracBsqrt2 A,bf a_1+fracsqrt2 A bf b_1Bright) +
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracB+2 A xsqrt2 A,bf a_1+fracsqrt2 A bf b_1B+2 A xright)+\
&&!!!!!!!!!!!!!!!!!!!!!!!!!!! -bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_1sqrt1+bf a_1^2,bf a_1+fracB(1+bf a_1^2)sqrt2 Abf b_1right) +
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_1sqrt1+bf a_1^2,bf a_1+frac(B+2 A x)(1+bf a_1^2)sqrt2 Abf b_1right)\
%
&&left.right)
endeqnarray
I can provide a code for testing the above expression if anyone is interested.
Now, in the particular case when all the diagonal elements of the matrix $bf A$ are equal unity and all the cross diagonal terms are equal to $rho$ where $0 le rho le 1$ then the result reads:
begineqnarray
&&P=\
&&frac2 pi ^3/2sqrt(1-rho )^3 (3 rho +1)
left( fracpi -3 arctanleft(sqrtfrac3 rho +1rho +1right)2 sqrtpi +6sqrtpi mathfrak J^(1,1)left( fracsqrtfrac32 rho sqrt(1-rho ) (3 rho +1),fracsqrt1-rho sqrt2 sqrt(1-rho ) (3 rho +1),sqrt3right)right)
endeqnarray
Below I plot the quantity $P$ as a function of $rho$. Note that the value $P(rho=0) = pi^2/4 simeq 2.4674$ as it is.
answered Jul 13 '17 at 16:19
PrzemoPrzemo
4,40811031
$endgroup$
add a comment |
$begingroup$
Thank you Przemo for your solution to the problem for $n=2, 3$. While I had no trouble following your derivation in 2D, I'm stuck with the derivation of your intermediate step for $n=3$. I mainly tried two approaches:
Completing the square in one variable, say $x$, leaves me with
$$int_mathbbR_+^2 mathrmdymathrmdz expleft(-frac12 fracmathrmdet,A_3mathrmdet,A_2z^2right) expleft(-frac12 fracmathrmdet, A_2a(y-m z)^2right) left[1 - mathrmerfleft(fraca_12y+a_13zsqrt2aright) right] $$
where $A_2=beginpmatrix a & a_12\ & bendpmatrix$, $A_3$ as you defined it, and $m$ is a function of the coefficients of the matrices. However, I do not know how to proceed from there: expanding the error function to do the integral in y, say, is a nightmare due to the constant term in z; I also did not find a way to do a coordinate transform à la $s=a_12y+a_13z$ or something similar.Indeed, your intermediate solution looks more like you were able to complete the square in two of the variables independently; but what happened to the cross-term? I cannot find a factorisation of the exponent that would allow me to complete two integrals over the half-line with only one variable left in the error function yielded by the integral.
Any help / hint would be greatly appreciated! Thank you in advance.
answered Dec 3 '18 at 18:09
workandheatworkandheat
587
$endgroup$
1
$begingroup$
@ workandheat, the key word here is orthant probability. If you search for this in connection with normal distribution, you get many results. But as I see even for cases with four variables people write long papers. So for more variables I doubt that there are exact results.
$endgroup$
– Karl
Dec 3 '18 at 18:23
$begingroup$
Thanks @Karl for the hint to "orthant" probability; coming from physics, I was not aware of this keyword. I appreciate this is a difficult problem in general, but for the moment I would just be happy to understand the derivation of Przemo's result...
$endgroup$
– workandheat
Dec 3 '18 at 18:31
$begingroup$
@workandheat Now the next step is to change variables $(y,z) rightarrow (u:=(a_1,2 y+a_1,3 z)/(sqrt2a),z)$. This changes the integration region from the first quadrant to an interior of an angle $0le u < infty$ and $0le z le cdots u$. Having done that we can still do the integral over $z$ by expressing it through error functions. Finally we end up with a one dimensional integral over $u$ from one Gaussian and a product of two error functions. We do the last integral using the result from the link provided.
$endgroup$
– Przemo
Jan 23 at 11:27
add a comment |
$begingroup$
Other names for this quantity are the "multivariate Gaussian cumulative distribution", the "normalization constant of the truncated normal distribution", "non-centered orthant probabilities", ...
There appears to be a rather extensive literature on this. See for example The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting and many citations therein, like this one
Here is a paper that has closed-form expressions for the orthant probabilities for $n=4$, under different sets of assumptions for the covariance matrix.
I will update this answer as I learn more about it
answered Feb 13 at 15:18
guillefixguillefix
295111
$endgroup$
add a comment |
$begingroup$
Here we provide an answer for $n=5$ in case when the underlying matrix $bf A$ has all diagonal entries equal to unity and all cross-diagonal entries equal to $rho$ where $0le rho le 1$. We have derived the result basically in the same way as in my previous answer above, meaning by firstly bringing the quadratic form to a square in one variable and integrating over that variable and then by successively integrating over the remaining variables and reducing the dimension of the integral.
Firstly let us note that the function $mathfrak J^(1,1)$ is defined as in my previous answer above and then
let us also define the following:
beginequation
mathfrak J^(2,1)left( (a_1,a_2), b, cright):= intlimits_0^infty frace^-1/2 xi^2sqrt2 pi cdot [prodlimits_j=1^2 erf(a_j xi)] cdot T(b xi, c) dxi
endequation
This function can be always reduced to di-logarithms as shown in An integral involving a Gaussian, error functions and the Owen's T function. .
Then we have:
begineqnarray
&&P=frac2^3/2 pisqrt(1-rho)^4(1+4 rho)left(right.\
&½sqrtpi mbox$left( 2 pi arctanleft(fracrho +1sqrt1-rho ^2right)-2 arctanleft(sqrtfrac4 rho +12 rho +1right) left(3 arctanleft(fracrho +1sqrt1-rho ^2right)+pi right)+arccos(rho )
left(pi -3 arctanleft(sqrtfrac4 rho +12 rho +1right)right)right)+$\
&& 6 pi^3/2 mathfrak J^(1,1)left( fracsqrtfrac32 rho sqrt4 rho +1,fracsqrtfracrho +14 rho +1sqrt2,sqrt3right)+\
&& 6 pi^3/2 mathfrak J^(2,1)left((fracsqrt4 rho +1sqrt6 rho ,fracrho +1sqrt2 sqrt1-rho ^2),fracsqrtrho +1sqrt6 rho ,sqrt3 right)+ \
&& 6 pi^3/2 mathfrak J^(2,1)left((sqrtfrac32,sqrtfrac4 rho +1rho +1),fracsqrt6 rho sqrtrho +1,fracrho +1sqrt1-rho ^2 right)\
left.right)
endeqnarray
Below I plot the quantity in question as a function of $rho$. Note that the value $P(rho=0)= (sqrtpi/sqrt2)^5 simeq 3.09243$ as it is.
answered Mar 14 at 17:27
PrzemoPrzemo
4,40811031
$endgroup$
add a comment |
Your Answer
StackExchange.ifUsing("editor", function ()
return StackExchange.using("mathjaxEditing", function ()
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix)
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
);
);
, "mathjax-editing");
StackExchange.ready(function()
var channelOptions =
tags: "".split(" "),
id: "69"
;
initTagRenderer("".split(" "), "".split(" "), channelOptions);
StackExchange.using("externalEditor", function()
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled)
StackExchange.using("snippets", function()
createEditor();
);
else
createEditor();
);
function createEditor()
StackExchange.prepareEditor(
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: true,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
imageUploader:
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
,
noCode: true, onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
);
);
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f869502%2fmultivariate-gaussian-integral-over-positive-reals%23new-answer', 'question_page');
);
Post as a guest
Required, but never shown
5 Answers
5
active
oldest
votes
5 Answers
5
active
oldest
votes
active
oldest
votes
active
oldest
votes
$begingroup$
The integral over (coordinate-wise) positive values appears in the treatment of dichotomized Gaussian distributions, so you might find the answer to your problem there. Relevant references would be:
- DR Cox, N Wermuth, Biometrika, 2002
- JH Macke, P Berens et al., Neural Computation, 2009
answered Jul 17 '14 at 14:51
flipsflips
515
$endgroup$
$begingroup$
Thanks! I looked into your referenced papers and it seems they approximate the integral in question under various assumptions, e.g. all covariances being identical.
$endgroup$
– le_m
Aug 4 '14 at 23:28
$begingroup$
As far as I can see, these only cite a solution for the case n=2, tracing back to Sheppard 1898
$endgroup$
– guillefix
Feb 13 at 15:09
add a comment |
$begingroup$
The integral over (coordinate-wise) positive values appears in the treatment of dichotomized Gaussian distributions, so you might find the answer to your problem there. Relevant references would be:
- DR Cox, N Wermuth, Biometrika, 2002
- JH Macke, P Berens et al., Neural Computation, 2009
answered Jul 17 '14 at 14:51
flipsflips
515
$endgroup$
$begingroup$
Thanks! I looked into your referenced papers and it seems they approximate the integral in question under various assumptions, e.g. all covariances being identical.
$endgroup$
– le_m
Aug 4 '14 at 23:28
$begingroup$
As far as I can see, these only cite a solution for the case n=2, tracing back to Sheppard 1898
$endgroup$
– guillefix
Feb 13 at 15:09
add a comment |
$begingroup$
The integral over (coordinate-wise) positive values appears in the treatment of dichotomized Gaussian distributions, so you might find the answer to your problem there. Relevant references would be:
- DR Cox, N Wermuth, Biometrika, 2002
- JH Macke, P Berens et al., Neural Computation, 2009
answered Jul 17 '14 at 14:51
flipsflips
515
$endgroup$
The integral over (coordinate-wise) positive values appears in the treatment of dichotomized Gaussian distributions, so you might find the answer to your problem there. Relevant references would be:
- DR Cox, N Wermuth, Biometrika, 2002
- JH Macke, P Berens et al., Neural Computation, 2009
answered Jul 17 '14 at 14:51
flipsflips
515
answered Jul 17 '14 at 14:51
flipsflips
515
answered Jul 17 '14 at 14:51
flipsflips
515
answered Jul 17 '14 at 14:51
flipsflips
515
515
$begingroup$
Thanks! I looked into your referenced papers and it seems they approximate the integral in question under various assumptions, e.g. all covariances being identical.
$endgroup$
– le_m
Aug 4 '14 at 23:28
$begingroup$
As far as I can see, these only cite a solution for the case n=2, tracing back to Sheppard 1898
$endgroup$
– guillefix
Feb 13 at 15:09
add a comment |
$begingroup$
Thanks! I looked into your referenced papers and it seems they approximate the integral in question under various assumptions, e.g. all covariances being identical.
$endgroup$
– le_m
Aug 4 '14 at 23:28
$begingroup$
As far as I can see, these only cite a solution for the case n=2, tracing back to Sheppard 1898
$endgroup$
– guillefix
Feb 13 at 15:09
$begingroup$
Thanks! I looked into your referenced papers and it seems they approximate the integral in question under various assumptions, e.g. all covariances being identical.
$endgroup$
– le_m
Aug 4 '14 at 23:28
$begingroup$
Thanks! I looked into your referenced papers and it seems they approximate the integral in question under various assumptions, e.g. all covariances being identical.
$endgroup$
– le_m
Aug 4 '14 at 23:28
$begingroup$
As far as I can see, these only cite a solution for the case n=2, tracing back to Sheppard 1898
$endgroup$
– guillefix
Feb 13 at 15:09
$begingroup$
As far as I can see, these only cite a solution for the case n=2, tracing back to Sheppard 1898
$endgroup$
– guillefix
Feb 13 at 15:09
add a comment |
$begingroup$
Let us compute the result in case $n=2$. Here the matrix reads $A=left(beginarrayrra & c\c& bendarrayright)$ .Therefore we have:
begineqnarray
P&=& intlimits_mathbb R_+^2 expleft-frac12left[sqrta(s_1+fracca s_2)right]^2 -frac12 fracb a-c^2a s_2^2right ds_1 ds_2\
&=&frac1sqrta sqrtfracpi2 intlimits_0^infty erfcleft(fraccsqrta fracs_2sqrt2 right)expleft-frac12(fracb a-c^2a)s_2^2 rightds_2\
&=&sqrtfracpi2 frac1sqrtb a-c^2 intlimits_0^infty erfc(fraccsqrtb a-c^2 fracs_2sqrt2) e^-frac12 s_2^2 ds_2\
&=& sqrtfracpi2 frac1sqrtb a-c^2 left(
sqrtfracpi2- sqrtfrac2pi arctan(fraccsqrtb a-c^2)right)\
&=& frac1sqrtb a-c^2 arctan(fracsqrtb a-c^2c)
endeqnarray
In the top line we completed the first integration variable to a square and in the second line we integrated over that variable. In the third line we changed variables accordingly . In the fourth line we integrated over the second variable by writing $erfc() = 1- erf()$ and then expanding the error function in a Taylor series and integrating term by term and finally in the last line we simplified the result.
Now, by doing similar calculations we obtained the following result in case $n=3$. Here $A=left(beginarrayrrra & a_12 & a_13\a_12& b&a_23\a_13&a_23&cendarrayright)$.
Firstly we have:
begineqnarray
&&vecs^(T).(A.vecs) = \
&&left(sqrta ( s_1 + fraca_1,2 s_2 + a_1,3 s_3a )right)^2 +
left( b- fraca_1,2^2aright) s_2^2 + left(c-fraca_1,3^2aright) s_3^2 + 2 left(a_2,3-fraca_1,2 a_1,3 aright) s_2 s_3
endeqnarray
Therefore integrating over $s_1$ gives:
begineqnarray
&&P=sqrtfracpi 2 frac1sqrta cdot \
&&intlimits_bf R^2 texterfcleft(fraca_1,2 s_2+a_1,3 s_3sqrt2 sqrtaright) cdot \
&&exp left[
-frac12 left(s_2^2 left(b-fraca_1,2^2aright)+2
s_2 s_3 left(a_2,3-fraca_1,2 a_1,3aright)+s_3^2 left(c-fraca_1,3^2aright)right)
right]
ds_2 ds_3=\
&&
fracsqrtpi a_1,2 intlimits_0^infty
texterfc(u) cdot expleft[-frac12u^2 (frac2 a ba_1,2^2 - 2)right]\
&&
intlimits_0^fracsqrt2 aa_1,3 u
exp left[-frac12 left(s_3 ufrac2 sqrt2 sqrta a_1,2 left(a_2,3-fracb
a_1,3a_1,2right)+
s_3^2fraca_1,3 a_1,2 left(fraca_1,3 ba_1,2+fraca_1,2 ca_1,3-2 a_2,3right)right)right] ds_3 du
endeqnarray
Now it is clear that we can do the integral over $s_3$ in the sense that we can express it through a difference of error functions.Denote $delta:=-2 a_1,2 a_1,3 a_2,3 +a_1,3^2 b +a_1,2^2 c$. Then we have
begineqnarray
&&P=fracpisqrt2sqrtdelta cdotintlimits_0^infty erfc(u) left( erfleft[fracsqrta(-a_1,3 a_2,3+a_1,2 c)a_1,3 sqrtdelta u right] - erfleft[ fracsqrta(a_1,2 a_2,3-a_1,3 b)a_1,2 sqrtdelta u right]right) e^-fracdet(A) delta u^2 du=\
&&fracpisqrt2 det(A)cdot \
&&
intlimits_0^infty erfcleft(u sqrtfracdeltadet(A)right)e^-u^2cdot \
&&left(-erfc(sqrta frac(-a_13a_23+a_12 c)a_13 sqrtdet(A) u)+erfc(sqrta frac(a_12a_23-a_13 b)a_12 sqrtdet(A) u)right) du \
&&=sqrtfracpi2 det(A)\
left[right.\
&&-arctanleft(fraca_13 sqrtdet(A)sqrta(-a_13a_23+a_12 c)right)+
arctanleft(fracsqrtc sqrtdet(A)-a_13 a_23 + a_12 cright)
\
&&+arctanleft(fraca_12 sqrtdet(A)sqrta (a_12 a_23 - a_13 b)right)-arctanleft(fracsqrtb sqrtdet(A)a_12 a_23 - a_13 bright)
left. right]\
&&=sqrtfracpi2 det(A)\
&&left[right.\
&&left.
arctanleft(frac(a_1,3-sqrta_1,1a_3,3)(a_1,3a_2,3-a_1,2a_3,3)sqrta_1,1 (a_1,3a_2,3-a_1,2a_3,3)^2+a_1,3 sqrta_3,3 det(A) sqrtdet(A)right)+right.\
&&left.
arctanleft(frac(a_1,2-sqrta_1,1a_2,2)(a_1,2a_2,3-a_1,3a_2,2)sqrta_1,1 (a_1,2a_2,3-a_1,3a_2,2)^2+a_1,2 sqrta_2,2 det(A) sqrtdet(A)right)
right]
endeqnarray
where in the last line we used An integral involving error functions and a Gaussian .
I also include a Mathematica code snippet that verifies all the steps involved:
(*3d*)
A =.; B =.; CC =.; A12 =.; A23 =.; A13 =.;
For[DDet = 0, True, ,
A, B, CC, A12, A23, A13 =
RandomReal[0, 1, 6, WorkingPrecision -> 50];
DDet = Det[A, A12, A13, A12, B, A23, A13, A23, CC];
If[DDet > 0, Break[]];
];
a = Sqrt[(-2 A12 A13 A23 + A13^2 B + A12^2 CC)/DDet];
b1, b2 = ( Sqrt[A] (-A13 A23 + A12 CC))/ Sqrt[DDet], (
Sqrt[A] (A12 A23 - A13 B))/ Sqrt[DDet];
AA1, AA2 = 2 Sqrt[2] Sqrt[
A] (( A23 A12 - A13 B)/A12^2), (-2 A12 A13 A23 + A13^2 B +
A12^2 CC)/A12^2;
DDet, a, b1, b2;
NIntegrate[
Exp[-1/2 (A s1^2 + B s2^2 + CC s3^2 + 2 A12 s1 s2 + 2 A23 s2 s3 +
2 A13 s1 s3)], s1, 0, Infinity, s2, 0, Infinity, s3, 0,
Infinity]
NIntegrate[
Exp[-1/2 ((Sqrt[A] (s1 + (A12 s2 + A13 s3)/A))^2 + (B -
A12^2/A) s2^2 + (CC - A13^2/A) s3^2 +
2 (A23 - A12 A13/A) s2 s3)], s1, 0, Infinity, s2, 0,
Infinity, s3, 0, Infinity]
NIntegrate[
1/Sqrt[A] Sqrt[
Pi/2] Erfc[(A12 s2 + A13 s3)/
Sqrt[2 A]] Exp[-1/
2 ((B - A12^2/A) s2^2 + (CC - A13^2/A) s3^2 +
2 (A23 - A12 A13/A) s2 s3)], s2, 0, Infinity, s3, 0,
Infinity]
Sqrt[Pi]/A12 NIntegrate[
Erfc[u] Exp[-1/
2 ( A13/A12 (-2 A23 + (A13 B)/A12 + CC A12/A13) s3^2 + (
2 Sqrt[2] Sqrt[A] )/
A12 ( A23 - ( A13 B)/A12) s3 u + (-2 + (2 A B)/
A12^2) u^2)], u, 0, Infinity, s3, 0, Sqrt[2 A]/A13 u]
Sqrt[Pi]/A12 NIntegrate[
Erfc[u] Exp[-1/2 (Sqrt[AA2] s3 + u/2 AA1/Sqrt[AA2])^2] Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity, s3, 0, Sqrt[2 A]/A13 u]
Sqrt[Pi]/(A12 Sqrt[AA2])
NIntegrate[
Erfc[u] Exp[-1/2 (s3)^2] Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity, s3,
u/2 AA1/Sqrt[AA2], ((A13 AA1 + 2 AA2 Sqrt[2] Sqrt[A]) u)/(
2 A13 Sqrt[AA2])]
Sqrt[Pi]/(A12 Sqrt[AA2]) Sqrt[[Pi]/2]
NIntegrate[
Erfc[u] (
Erf[(A13 AA1 + 2 AA2 Sqrt[2] Sqrt[A])/(2 A13 Sqrt[2] Sqrt[AA2])
u] - Erf[AA1/(2 Sqrt[2] Sqrt[AA2]) u]) Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity]
Pi/Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC] Sqrt[1/2]
NIntegrate[
Erfc[u] (
Erf[( Sqrt[A] (-A13 A23 + A12 CC) u)/(
A13 Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC])] -
Erf[(Sqrt[A] (A12 A23 - A13 B) u)/(
A12 Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC])]) Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity]
Pi/ Sqrt[-2 A12 A13 A23 + A13^2 B +
A12^2 CC] Sqrt[1/2] a NIntegrate[
Erfc[a u] (
Erf[( Sqrt[A] (-A13 A23 + A12 CC) u)/(A13 Sqrt[DDet])] -
Erf[(Sqrt[A] (A12 A23 - A13 B) u)/(A12 Sqrt[DDet])]) Exp[-
u^2], u, 0, Infinity]
Pi/Sqrt[2 DDet] NIntegrate[(Erfc[u a]) Exp[-u^2] (Erf[b1/A13 u] -
Erf[b2/A12 u]), u, 0, Infinity]
Sqrt[Pi]/Sqrt[
2 DDet] (ArcTan[ Sqrt[A]/A13 (-A13 A23 + A12 CC)/ Sqrt[DDet]] -
ArcTan[1/ Sqrt[CC] (-A13 A23 + A12 CC)/ Sqrt[DDet]] -
ArcTan[ Sqrt[A]/A12 (A12 A23 - A13 B)/ Sqrt[DDet]] +
ArcTan[ 1/Sqrt[B] (A12 A23 - A13 B)/ Sqrt[DDet]])
-(Sqrt[Pi]/
Sqrt[2 DDet]) (ArcTan[(A13 Sqrt[DDet])/(
Sqrt[A] (-A13 A23 + A12 CC))] -
ArcTan[(Sqrt[CC] Sqrt[DDet])/(-A13 A23 + A12 CC)] -
ArcTan[(A12 Sqrt[DDet])/(Sqrt[A] (A12 A23 - A13 B))] +
ArcTan[(Sqrt[B] Sqrt[DDet])/(A12 A23 - A13 B)])
Sqrt[Pi]/Sqrt[
2 DDet] (ArcTan[((A13 - Sqrt[A] Sqrt[CC]) (A13 A23 - A12 CC) Sqrt[
DDet])/(Sqrt[A] (A13 A23 - A12 CC)^2 + A13 Sqrt[CC] DDet)] +
ArcTan[((A12 - Sqrt[A] Sqrt[B]) (A12 A23 - A13 B) Sqrt[DDet])/(
Sqrt[A] (A12 A23 - A13 B)^2 + A12 Sqrt[B] DDet)])
Update: Now let us take a look at the $n=4$ case. In here:
beginequation
bf A=left(
beginarrayrrrr
a & a_1,2 & a_1,3 & a_1,4 \
a_1,2 & b & a_2,3 & a_2,4 \
a_1,3 & a_2,3 & c & a_3,4 \
a_1,4 & a_2,4 & a_3,4 & d
endarray
right)
endequation
then by doing basically the same computations as a above we managed to reduce the integral in question to a following two dimensional integral. We have:
begineqnarray
&&P= \
&&!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!fracpisqrt2 delta intlimits_0^infty intlimits_0^fracsqrt2 aa_1,2 u
erfc[u] cdot expleft[fracmathfrak A_0,0 u^2 + mathfrak A_1,0 u s_2 +mathfrak A_1,1 s_2^22 deltaright] cdot left( erf[fracmathfrak B_1 u + mathfrak B_2 s_2a_1,3 sqrt2 delta] + erf[fracmathfrak C_1 u + mathfrak C_2 s_2a_1,4 sqrt2 delta]right) d s_2 du=\
&&frac2 imath pi^3/2sqrtmathfrak A_1,1
intlimits_0^infty erfc[u] expfrac4 mathfrak A_0,0 mathfrak A_1,1 - mathfrak A_1,0^28delta mathfrak A_1,1 u^2cdot\
&&
left[right.\
&&left.
left.Tleft(frac(mathfrak A_1,0+xi) u2imath sqrtmathfrak A_1,1 delta,
fracimath mathfrak B_2a_1,3 sqrtmathfrak A_1,1,fracu(2mathfrak A_1,1 mathfrak B_1 - mathfrak A_1,0 mathfrak B_2)2sqrtdelta a_1,3 mathfrak A_1,1right)right|_frac2mathfrak A_1,1 sqrt2 aa_1,2^0 +right.\
&&left.
left.Tleft(frac(mathfrak A_1,0+xi) u2imath sqrtmathfrak A_1,1 delta,
fracimath mathfrak C_2a_1,3 sqrtmathfrak A_1,1,fracu(2mathfrak A_1,1 mathfrak C_1 - mathfrak A_1,0 mathfrak C_2)2sqrtdelta a_1,3 mathfrak A_1,1right)right|_frac2mathfrak A_1,1 sqrt2 aa_1,2^0 +right.\
&&left.
right] du quad (i)
endeqnarray
where $T(cdot,cdot,cdot)$ is the generalized Owen's T function Generalized Owen's T function and
begineqnarray
delta&:=&a_1,3(a_1,3 d-a_1,4 a_3,4) + a_1,4(a_1,4 c- a_1,3 a_3,4)\
mathfrak A_0,0&:=&2 a left(a_3,4^2-c dright)+2 a_1,4 (a_1,4 c-a_1,3 a_3,4)+2 a_1,3 (a_1,3 d-a_1,4 a_3,4)\
mathfrak A_1,0&:=&2 sqrt2 sqrta left(a_1,2 left(c d-a_3,4^2right)+a_1,3
(a_2,4 a_3,4-a_2,3 d)+a_1,4 (a_2,3 a_3,4-a_2,4 c)right)\
mathfrak A_1,1&:=&a_1,2^2 left(a_3,4^2-c dright)+2 a_1,2 a_1,3 (a_2,3 d-a_2,4 a_3,4)+2 a_1,2 a_1,4
(a_2,4 c-a_2,3 a_3,4)+a_1,3^2 left(a_2,4^2-b dright)+2 a_1,3 a_1,4 (a_3,4 b-a_2,3 a_2,4)+a_1,4^2 left(a_2,3^2-b cright)\
hline\
mathfrak B_1&:=&sqrt2 sqrta (a_1,4
c-a_1,3 a_3,4)\
mathfrak B_2&:=&a_1,2 (a_1,3 a_3,4-a_1,4 c)+a_1,3 (a_1,4 a_2,3-a_1,3 a_2,4)\
mathfrak C_1&:=&sqrt2 sqrta (a_1,3 d-a_1,4 a_3,4)\
mathfrak C_2&:=&a_1,2 (a_1,4
a_3,4-a_1,3 d)+a_1,4 (a_1,3 a_2,4-a_1,4 a_2,3)
endeqnarray
nu = 4; Clear[T]; Clear[a]; x =.;
(*a0.dat, a1.dat or a2.dat*)
mat = << "a0.dat";
a, b, c, d, a12, a13, a14, a23, a24, a34 = mat[[1, 1]],
mat[[2, 2]], mat[[3, 3]], mat[[4, 4]], mat[[1, 2]], mat[[1, 3]],
mat[[1, 4]], mat[[2, 3]], mat[[2, 4]], mat[[3, 4]];
dd, A00, A10,
A11 = -2 a13 a14 a34 + a14^2 c + a13^2 d, -4 a13 a14 a34 +
2 a a34^2 + 2 a14^2 c + 2 a13^2 d - 2 a c d,
2 Sqrt[2] Sqrt[a] a14 a23 a34 + 2 Sqrt[2] Sqrt[a] a13 a24 a34 -
2 Sqrt[2] Sqrt[a] a12 a34^2 - 2 Sqrt[2] Sqrt[a] a14 a24 c -
2 Sqrt[2] Sqrt[a] a13 a23 d + 2 Sqrt[2] Sqrt[a] a12 c d,
a14^2 a23^2 - 2 a13 a14 a23 a24 + a13^2 a24^2 -
2 a12 a14 a23 a34 - 2 a12 a13 a24 a34 + a12^2 a34^2 +
2 a13 a14 a34 b + 2 a12 a14 a24 c - a14^2 b c + 2 a12 a13 a23 d -
a13^2 b d - a12^2 c d;
B1, B2, C1,
C2 = Sqrt[2] Sqrt[
a] (-a13 a34 + a14 c), (a13 a14 a23 - a13^2 a24 + a12 a13 a34 -
a12 a14 c),
Sqrt[2] Sqrt[
a] (-a14 a34 + a13 d), (-a14^2 a23 + a13 a14 a24 + a12 a14 a34 -
a12 a13 d);
NIntegrate[
Exp[-1/2 Sum[mat[[i, j]] s[i] s[j], i, 1, nu, j, 1, nu]],
Evaluate[Sequence @@ Table[s[eta], 0, Infinity, eta, 1, nu]]]
Sqrt[[Pi]/(2 a)]
NIntegrate[
Erfc[(a12 s[2] + a13 s[3] + a14 s[4])/Sqrt[
2 a]] Exp[-1/
2 ((-(a12^2/a) + b) s[2]^2 + (-(a13^2/a) + c) s[
3]^2 + (-(a14^2/a) + d) s[4]^2 +
2 (-(( a13 a14)/a) + a34) s[3] s[4] +
2 (-(( a12 a13)/a) + a23) s[2] s[3] +
2 (-(( a12 a14)/a) + a24) s[2] s[4])],
Evaluate[Sequence @@ Table[s[eta], 0, Infinity, eta, 2, nu]]]
Sqrt[[Pi]]
1/a14 NIntegrate[
Erfc[u] Exp[(
2 a14 a24 s[2] (-Sqrt[2] Sqrt[a] u + a12 s[2]) -
d (2 a u^2 - 2 Sqrt[2] Sqrt[a] a12 u s[2] + a12^2 s[2]^2) +
a14^2 (2 u^2 - b s[2]^2))/(
2 a14^2) + ((Sqrt[2] Sqrt[
a] (-a14 a34 + a13 d) u + (-a14^2 a23 + a13 a14 a24 +
a12 a14 a34 - a12 a13 d) s[2]) s[3])/
a14^2 - ((-2 a13 a14 a34 + a14^2 c + a13^2 d) s[3]^2)/(
2 a14^2)], u, 0, Infinity, s[2], 0,
Sqrt[2] Sqrt[a]/a12 u, s[3], 0, (Sqrt[2 a] u - a12 s[2])/a13]
Pi/Sqrt[2 dd]
NIntegrate[
Erfc[u] Exp[(A00 u^2 + A10 u s[2] + A11 s[2]^2)/(
2 (dd))] (Erf[(B1 u + B2 s[2])/( a13 Sqrt[2 dd])] +
Erf[(C1 u + C2 s[2])/( a14^1 Sqrt[2 dd])]), u, 0,
Infinity, s[2], 0, Sqrt[2] Sqrt[a]/a12 u]
Now, I will provide the result. Note that the only assumptions on the underlying matrix $bf A$ are that it is symmetric and that its elements are non-negative. Firstly let us define:
begineqnarray
&&mathfrak J^(1,1)(a,b,c)= frac1pi^2cdot left(right.\
&&left.
-frac18
sumlimits_i=1^4 sumlimits_j=1^4
(-1)^j-1+lfloor fraci-12 rfloor
%
mathfrak F^(1,fracsqrt1+2 a^2+b^2 - sqrt2 asqrt1+b^2)_fraci sqrtb^2 c^2+b^2+1 (-1)^leftlfloor
fracj-12rightrfloor +i b c (-1)^jsqrtb^2+1,-fracb
(-1)^i+i (-1)^leftlceil fraci-12rightrceil
sqrtb^2+1
%
right. \
&&left.
right)quad (ii)
endeqnarray
where $mathfrak F^(A,B)_a,b$ is related to di-logarithms and is defined in An integral involving a Gaussian, error functions and the Owen's T function. .
Then we define another function as follows:
beginequation
bar mathfrak J^(1,1)(a,b,c):= fracpi2 arctanleft[ fracsqrt2 a csqrt2 a+b^2(1+c^2)right] - fracpi2 arctanleft[ cright] - 2 pi^2 mathfrak J^(1,1)(frac1sqrt2 a,fracbsqrt2 a,c)
endequation
and then the following quantities that depend on the underlying matrix. We have:
begineqnarray
delta&:=& a_3,3 a_4,1^2 - 2 a_3,1 a_3,4 a_4,1 + a_4,4 a_3,1^2\
W&:=&left(a_3,3 a_4,4-a_3,4^2right) a_1,2^2+2 a_1,4 (a_2,3
a_3,4-a_2,4 a_3,3) a_1,2+2 a_1,3 (a_2,4 a_3,4-a_2,3 a_4,4)
a_1,2+a_1,4^2 left(a_2,2 a_3,3-a_2,3^2right)+2 a_1,3 a_1,4
(a_2,3 a_2,4-a_2,2 a_3,4)+a_1,3^2 left(a_2,2
a_4,4-a_2,4^2right)\
W_1&:=&2 sqrta_1,1 left(a_1,4 (a_2,4 a_3,3-a_2,3 a_3,4)+a_1,3 (a_2,3
a_4,4-a_2,4 a_3,4)+a_1,2 left(a_3,4^2-a_3,3 a_4,4right)right)\
%
v_1&:=&frac1a_4,1 sqrtdelta left( sqrta_1,1(a_3,4 a_4,1 - a_3,1 a_4,4),-a_2,4 a_3,1 a_4,1 + a_2,3 a_4,1^2+a_2,1(-a_3,4 a_4,1+a_3,1a_4,4)right)\
v_2&:=&-frac1a_3,1 sqrtdelta left(sqrta_1,1(a_3,4 a_3,1 - a_4,1 a_3,3),-a_3,1 a_3,2 a_4,1 +a_2,4 a_3,1^2 + a_2,1(-a_3,4 a_3,1+a_4,1a_3,3) right)\
%
left( A, B right)&:=& frac1delta left( W,W_1 right)\
left( bf a_1,bf a_2 right)&:=& frac1sqrtA left(v_1(2),v_2(2) right)\
bf b_1&:=& sqrt2 v_1(1) - fracBsqrt2 A v_1(2)\
bf b_2&:=& sqrt2 v_2(1) - fracBsqrt2 A v_2(2)\
x&:=& fracsqrta_1,1a_2,1
endeqnarray
Then the result reads:
begineqnarray
&&P=frac1det(bf A) left(right.\
%
&& bar mathfrak J^(1,1)left( fracdet(bf A)W,fracBsqrt2 A,bf a_2+fracsqrt2 A bf b_2Bright) -
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracB+2 A xsqrt2 A,bf a_2+fracsqrt2 A bf b_2B+2 A xright)+\
&&!!!!!!!!!!!!!!!!!!!!!!!!!!! bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_2sqrt1+bf a_2^2,bf a_2+fracB(1+bf a_2^2)sqrt2 Abf b_2right) -
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_2sqrt1+bf a_2^2,bf a_2+frac(B+2 A x)(1+bf a_2^2)sqrt2 Abf b_2right)+\
%
&& -bar mathfrak J^(1,1)left( fracdet(bf A)W,fracBsqrt2 A,bf a_1+fracsqrt2 A bf b_1Bright) +
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracB+2 A xsqrt2 A,bf a_1+fracsqrt2 A bf b_1B+2 A xright)+\
&&!!!!!!!!!!!!!!!!!!!!!!!!!!! -bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_1sqrt1+bf a_1^2,bf a_1+fracB(1+bf a_1^2)sqrt2 Abf b_1right) +
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_1sqrt1+bf a_1^2,bf a_1+frac(B+2 A x)(1+bf a_1^2)sqrt2 Abf b_1right)\
%
&&left.right)
endeqnarray
I can provide a code for testing the above expression if anyone is interested.
Now, in the particular case when all the diagonal elements of the matrix $bf A$ are equal unity and all the cross diagonal terms are equal to $rho$ where $0 le rho le 1$ then the result reads:
begineqnarray
&&P=\
&&frac2 pi ^3/2sqrt(1-rho )^3 (3 rho +1)
left( fracpi -3 arctanleft(sqrtfrac3 rho +1rho +1right)2 sqrtpi +6sqrtpi mathfrak J^(1,1)left( fracsqrtfrac32 rho sqrt(1-rho ) (3 rho +1),fracsqrt1-rho sqrt2 sqrt(1-rho ) (3 rho +1),sqrt3right)right)
endeqnarray
Below I plot the quantity $P$ as a function of $rho$. Note that the value $P(rho=0) = pi^2/4 simeq 2.4674$ as it is.
answered Jul 13 '17 at 16:19
PrzemoPrzemo
4,40811031
$endgroup$
add a comment |
$begingroup$
Let us compute the result in case $n=2$. Here the matrix reads $A=left(beginarrayrra & c\c& bendarrayright)$ .Therefore we have:
begineqnarray
P&=& intlimits_mathbb R_+^2 expleft-frac12left[sqrta(s_1+fracca s_2)right]^2 -frac12 fracb a-c^2a s_2^2right ds_1 ds_2\
&=&frac1sqrta sqrtfracpi2 intlimits_0^infty erfcleft(fraccsqrta fracs_2sqrt2 right)expleft-frac12(fracb a-c^2a)s_2^2 rightds_2\
&=&sqrtfracpi2 frac1sqrtb a-c^2 intlimits_0^infty erfc(fraccsqrtb a-c^2 fracs_2sqrt2) e^-frac12 s_2^2 ds_2\
&=& sqrtfracpi2 frac1sqrtb a-c^2 left(
sqrtfracpi2- sqrtfrac2pi arctan(fraccsqrtb a-c^2)right)\
&=& frac1sqrtb a-c^2 arctan(fracsqrtb a-c^2c)
endeqnarray
In the top line we completed the first integration variable to a square and in the second line we integrated over that variable. In the third line we changed variables accordingly . In the fourth line we integrated over the second variable by writing $erfc() = 1- erf()$ and then expanding the error function in a Taylor series and integrating term by term and finally in the last line we simplified the result.
Now, by doing similar calculations we obtained the following result in case $n=3$. Here $A=left(beginarrayrrra & a_12 & a_13\a_12& b&a_23\a_13&a_23&cendarrayright)$.
Firstly we have:
begineqnarray
&&vecs^(T).(A.vecs) = \
&&left(sqrta ( s_1 + fraca_1,2 s_2 + a_1,3 s_3a )right)^2 +
left( b- fraca_1,2^2aright) s_2^2 + left(c-fraca_1,3^2aright) s_3^2 + 2 left(a_2,3-fraca_1,2 a_1,3 aright) s_2 s_3
endeqnarray
Therefore integrating over $s_1$ gives:
begineqnarray
&&P=sqrtfracpi 2 frac1sqrta cdot \
&&intlimits_bf R^2 texterfcleft(fraca_1,2 s_2+a_1,3 s_3sqrt2 sqrtaright) cdot \
&&exp left[
-frac12 left(s_2^2 left(b-fraca_1,2^2aright)+2
s_2 s_3 left(a_2,3-fraca_1,2 a_1,3aright)+s_3^2 left(c-fraca_1,3^2aright)right)
right]
ds_2 ds_3=\
&&
fracsqrtpi a_1,2 intlimits_0^infty
texterfc(u) cdot expleft[-frac12u^2 (frac2 a ba_1,2^2 - 2)right]\
&&
intlimits_0^fracsqrt2 aa_1,3 u
exp left[-frac12 left(s_3 ufrac2 sqrt2 sqrta a_1,2 left(a_2,3-fracb
a_1,3a_1,2right)+
s_3^2fraca_1,3 a_1,2 left(fraca_1,3 ba_1,2+fraca_1,2 ca_1,3-2 a_2,3right)right)right] ds_3 du
endeqnarray
Now it is clear that we can do the integral over $s_3$ in the sense that we can express it through a difference of error functions.Denote $delta:=-2 a_1,2 a_1,3 a_2,3 +a_1,3^2 b +a_1,2^2 c$. Then we have
begineqnarray
&&P=fracpisqrt2sqrtdelta cdotintlimits_0^infty erfc(u) left( erfleft[fracsqrta(-a_1,3 a_2,3+a_1,2 c)a_1,3 sqrtdelta u right] - erfleft[ fracsqrta(a_1,2 a_2,3-a_1,3 b)a_1,2 sqrtdelta u right]right) e^-fracdet(A) delta u^2 du=\
&&fracpisqrt2 det(A)cdot \
&&
intlimits_0^infty erfcleft(u sqrtfracdeltadet(A)right)e^-u^2cdot \
&&left(-erfc(sqrta frac(-a_13a_23+a_12 c)a_13 sqrtdet(A) u)+erfc(sqrta frac(a_12a_23-a_13 b)a_12 sqrtdet(A) u)right) du \
&&=sqrtfracpi2 det(A)\
left[right.\
&&-arctanleft(fraca_13 sqrtdet(A)sqrta(-a_13a_23+a_12 c)right)+
arctanleft(fracsqrtc sqrtdet(A)-a_13 a_23 + a_12 cright)
\
&&+arctanleft(fraca_12 sqrtdet(A)sqrta (a_12 a_23 - a_13 b)right)-arctanleft(fracsqrtb sqrtdet(A)a_12 a_23 - a_13 bright)
left. right]\
&&=sqrtfracpi2 det(A)\
&&left[right.\
&&left.
arctanleft(frac(a_1,3-sqrta_1,1a_3,3)(a_1,3a_2,3-a_1,2a_3,3)sqrta_1,1 (a_1,3a_2,3-a_1,2a_3,3)^2+a_1,3 sqrta_3,3 det(A) sqrtdet(A)right)+right.\
&&left.
arctanleft(frac(a_1,2-sqrta_1,1a_2,2)(a_1,2a_2,3-a_1,3a_2,2)sqrta_1,1 (a_1,2a_2,3-a_1,3a_2,2)^2+a_1,2 sqrta_2,2 det(A) sqrtdet(A)right)
right]
endeqnarray
where in the last line we used An integral involving error functions and a Gaussian .
I also include a Mathematica code snippet that verifies all the steps involved:
(*3d*)
A =.; B =.; CC =.; A12 =.; A23 =.; A13 =.;
For[DDet = 0, True, ,
A, B, CC, A12, A23, A13 =
RandomReal[0, 1, 6, WorkingPrecision -> 50];
DDet = Det[A, A12, A13, A12, B, A23, A13, A23, CC];
If[DDet > 0, Break[]];
];
a = Sqrt[(-2 A12 A13 A23 + A13^2 B + A12^2 CC)/DDet];
b1, b2 = ( Sqrt[A] (-A13 A23 + A12 CC))/ Sqrt[DDet], (
Sqrt[A] (A12 A23 - A13 B))/ Sqrt[DDet];
AA1, AA2 = 2 Sqrt[2] Sqrt[
A] (( A23 A12 - A13 B)/A12^2), (-2 A12 A13 A23 + A13^2 B +
A12^2 CC)/A12^2;
DDet, a, b1, b2;
NIntegrate[
Exp[-1/2 (A s1^2 + B s2^2 + CC s3^2 + 2 A12 s1 s2 + 2 A23 s2 s3 +
2 A13 s1 s3)], s1, 0, Infinity, s2, 0, Infinity, s3, 0,
Infinity]
NIntegrate[
Exp[-1/2 ((Sqrt[A] (s1 + (A12 s2 + A13 s3)/A))^2 + (B -
A12^2/A) s2^2 + (CC - A13^2/A) s3^2 +
2 (A23 - A12 A13/A) s2 s3)], s1, 0, Infinity, s2, 0,
Infinity, s3, 0, Infinity]
NIntegrate[
1/Sqrt[A] Sqrt[
Pi/2] Erfc[(A12 s2 + A13 s3)/
Sqrt[2 A]] Exp[-1/
2 ((B - A12^2/A) s2^2 + (CC - A13^2/A) s3^2 +
2 (A23 - A12 A13/A) s2 s3)], s2, 0, Infinity, s3, 0,
Infinity]
Sqrt[Pi]/A12 NIntegrate[
Erfc[u] Exp[-1/
2 ( A13/A12 (-2 A23 + (A13 B)/A12 + CC A12/A13) s3^2 + (
2 Sqrt[2] Sqrt[A] )/
A12 ( A23 - ( A13 B)/A12) s3 u + (-2 + (2 A B)/
A12^2) u^2)], u, 0, Infinity, s3, 0, Sqrt[2 A]/A13 u]
Sqrt[Pi]/A12 NIntegrate[
Erfc[u] Exp[-1/2 (Sqrt[AA2] s3 + u/2 AA1/Sqrt[AA2])^2] Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity, s3, 0, Sqrt[2 A]/A13 u]
Sqrt[Pi]/(A12 Sqrt[AA2])
NIntegrate[
Erfc[u] Exp[-1/2 (s3)^2] Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity, s3,
u/2 AA1/Sqrt[AA2], ((A13 AA1 + 2 AA2 Sqrt[2] Sqrt[A]) u)/(
2 A13 Sqrt[AA2])]
Sqrt[Pi]/(A12 Sqrt[AA2]) Sqrt[[Pi]/2]
NIntegrate[
Erfc[u] (
Erf[(A13 AA1 + 2 AA2 Sqrt[2] Sqrt[A])/(2 A13 Sqrt[2] Sqrt[AA2])
u] - Erf[AA1/(2 Sqrt[2] Sqrt[AA2]) u]) Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity]
Pi/Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC] Sqrt[1/2]
NIntegrate[
Erfc[u] (
Erf[( Sqrt[A] (-A13 A23 + A12 CC) u)/(
A13 Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC])] -
Erf[(Sqrt[A] (A12 A23 - A13 B) u)/(
A12 Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC])]) Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity]
Pi/ Sqrt[-2 A12 A13 A23 + A13^2 B +
A12^2 CC] Sqrt[1/2] a NIntegrate[
Erfc[a u] (
Erf[( Sqrt[A] (-A13 A23 + A12 CC) u)/(A13 Sqrt[DDet])] -
Erf[(Sqrt[A] (A12 A23 - A13 B) u)/(A12 Sqrt[DDet])]) Exp[-
u^2], u, 0, Infinity]
Pi/Sqrt[2 DDet] NIntegrate[(Erfc[u a]) Exp[-u^2] (Erf[b1/A13 u] -
Erf[b2/A12 u]), u, 0, Infinity]
Sqrt[Pi]/Sqrt[
2 DDet] (ArcTan[ Sqrt[A]/A13 (-A13 A23 + A12 CC)/ Sqrt[DDet]] -
ArcTan[1/ Sqrt[CC] (-A13 A23 + A12 CC)/ Sqrt[DDet]] -
ArcTan[ Sqrt[A]/A12 (A12 A23 - A13 B)/ Sqrt[DDet]] +
ArcTan[ 1/Sqrt[B] (A12 A23 - A13 B)/ Sqrt[DDet]])
-(Sqrt[Pi]/
Sqrt[2 DDet]) (ArcTan[(A13 Sqrt[DDet])/(
Sqrt[A] (-A13 A23 + A12 CC))] -
ArcTan[(Sqrt[CC] Sqrt[DDet])/(-A13 A23 + A12 CC)] -
ArcTan[(A12 Sqrt[DDet])/(Sqrt[A] (A12 A23 - A13 B))] +
ArcTan[(Sqrt[B] Sqrt[DDet])/(A12 A23 - A13 B)])
Sqrt[Pi]/Sqrt[
2 DDet] (ArcTan[((A13 - Sqrt[A] Sqrt[CC]) (A13 A23 - A12 CC) Sqrt[
DDet])/(Sqrt[A] (A13 A23 - A12 CC)^2 + A13 Sqrt[CC] DDet)] +
ArcTan[((A12 - Sqrt[A] Sqrt[B]) (A12 A23 - A13 B) Sqrt[DDet])/(
Sqrt[A] (A12 A23 - A13 B)^2 + A12 Sqrt[B] DDet)])
Update: Now let us take a look at the $n=4$ case. In here:
beginequation
bf A=left(
beginarrayrrrr
a & a_1,2 & a_1,3 & a_1,4 \
a_1,2 & b & a_2,3 & a_2,4 \
a_1,3 & a_2,3 & c & a_3,4 \
a_1,4 & a_2,4 & a_3,4 & d
endarray
right)
endequation
then by doing basically the same computations as a above we managed to reduce the integral in question to a following two dimensional integral. We have:
begineqnarray
&&P= \
&&!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!fracpisqrt2 delta intlimits_0^infty intlimits_0^fracsqrt2 aa_1,2 u
erfc[u] cdot expleft[fracmathfrak A_0,0 u^2 + mathfrak A_1,0 u s_2 +mathfrak A_1,1 s_2^22 deltaright] cdot left( erf[fracmathfrak B_1 u + mathfrak B_2 s_2a_1,3 sqrt2 delta] + erf[fracmathfrak C_1 u + mathfrak C_2 s_2a_1,4 sqrt2 delta]right) d s_2 du=\
&&frac2 imath pi^3/2sqrtmathfrak A_1,1
intlimits_0^infty erfc[u] expfrac4 mathfrak A_0,0 mathfrak A_1,1 - mathfrak A_1,0^28delta mathfrak A_1,1 u^2cdot\
&&
left[right.\
&&left.
left.Tleft(frac(mathfrak A_1,0+xi) u2imath sqrtmathfrak A_1,1 delta,
fracimath mathfrak B_2a_1,3 sqrtmathfrak A_1,1,fracu(2mathfrak A_1,1 mathfrak B_1 - mathfrak A_1,0 mathfrak B_2)2sqrtdelta a_1,3 mathfrak A_1,1right)right|_frac2mathfrak A_1,1 sqrt2 aa_1,2^0 +right.\
&&left.
left.Tleft(frac(mathfrak A_1,0+xi) u2imath sqrtmathfrak A_1,1 delta,
fracimath mathfrak C_2a_1,3 sqrtmathfrak A_1,1,fracu(2mathfrak A_1,1 mathfrak C_1 - mathfrak A_1,0 mathfrak C_2)2sqrtdelta a_1,3 mathfrak A_1,1right)right|_frac2mathfrak A_1,1 sqrt2 aa_1,2^0 +right.\
&&left.
right] du quad (i)
endeqnarray
where $T(cdot,cdot,cdot)$ is the generalized Owen's T function Generalized Owen's T function and
begineqnarray
delta&:=&a_1,3(a_1,3 d-a_1,4 a_3,4) + a_1,4(a_1,4 c- a_1,3 a_3,4)\
mathfrak A_0,0&:=&2 a left(a_3,4^2-c dright)+2 a_1,4 (a_1,4 c-a_1,3 a_3,4)+2 a_1,3 (a_1,3 d-a_1,4 a_3,4)\
mathfrak A_1,0&:=&2 sqrt2 sqrta left(a_1,2 left(c d-a_3,4^2right)+a_1,3
(a_2,4 a_3,4-a_2,3 d)+a_1,4 (a_2,3 a_3,4-a_2,4 c)right)\
mathfrak A_1,1&:=&a_1,2^2 left(a_3,4^2-c dright)+2 a_1,2 a_1,3 (a_2,3 d-a_2,4 a_3,4)+2 a_1,2 a_1,4
(a_2,4 c-a_2,3 a_3,4)+a_1,3^2 left(a_2,4^2-b dright)+2 a_1,3 a_1,4 (a_3,4 b-a_2,3 a_2,4)+a_1,4^2 left(a_2,3^2-b cright)\
hline\
mathfrak B_1&:=&sqrt2 sqrta (a_1,4
c-a_1,3 a_3,4)\
mathfrak B_2&:=&a_1,2 (a_1,3 a_3,4-a_1,4 c)+a_1,3 (a_1,4 a_2,3-a_1,3 a_2,4)\
mathfrak C_1&:=&sqrt2 sqrta (a_1,3 d-a_1,4 a_3,4)\
mathfrak C_2&:=&a_1,2 (a_1,4
a_3,4-a_1,3 d)+a_1,4 (a_1,3 a_2,4-a_1,4 a_2,3)
endeqnarray
nu = 4; Clear[T]; Clear[a]; x =.;
(*a0.dat, a1.dat or a2.dat*)
mat = << "a0.dat";
a, b, c, d, a12, a13, a14, a23, a24, a34 = mat[[1, 1]],
mat[[2, 2]], mat[[3, 3]], mat[[4, 4]], mat[[1, 2]], mat[[1, 3]],
mat[[1, 4]], mat[[2, 3]], mat[[2, 4]], mat[[3, 4]];
dd, A00, A10,
A11 = -2 a13 a14 a34 + a14^2 c + a13^2 d, -4 a13 a14 a34 +
2 a a34^2 + 2 a14^2 c + 2 a13^2 d - 2 a c d,
2 Sqrt[2] Sqrt[a] a14 a23 a34 + 2 Sqrt[2] Sqrt[a] a13 a24 a34 -
2 Sqrt[2] Sqrt[a] a12 a34^2 - 2 Sqrt[2] Sqrt[a] a14 a24 c -
2 Sqrt[2] Sqrt[a] a13 a23 d + 2 Sqrt[2] Sqrt[a] a12 c d,
a14^2 a23^2 - 2 a13 a14 a23 a24 + a13^2 a24^2 -
2 a12 a14 a23 a34 - 2 a12 a13 a24 a34 + a12^2 a34^2 +
2 a13 a14 a34 b + 2 a12 a14 a24 c - a14^2 b c + 2 a12 a13 a23 d -
a13^2 b d - a12^2 c d;
B1, B2, C1,
C2 = Sqrt[2] Sqrt[
a] (-a13 a34 + a14 c), (a13 a14 a23 - a13^2 a24 + a12 a13 a34 -
a12 a14 c),
Sqrt[2] Sqrt[
a] (-a14 a34 + a13 d), (-a14^2 a23 + a13 a14 a24 + a12 a14 a34 -
a12 a13 d);
NIntegrate[
Exp[-1/2 Sum[mat[[i, j]] s[i] s[j], i, 1, nu, j, 1, nu]],
Evaluate[Sequence @@ Table[s[eta], 0, Infinity, eta, 1, nu]]]
Sqrt[[Pi]/(2 a)]
NIntegrate[
Erfc[(a12 s[2] + a13 s[3] + a14 s[4])/Sqrt[
2 a]] Exp[-1/
2 ((-(a12^2/a) + b) s[2]^2 + (-(a13^2/a) + c) s[
3]^2 + (-(a14^2/a) + d) s[4]^2 +
2 (-(( a13 a14)/a) + a34) s[3] s[4] +
2 (-(( a12 a13)/a) + a23) s[2] s[3] +
2 (-(( a12 a14)/a) + a24) s[2] s[4])],
Evaluate[Sequence @@ Table[s[eta], 0, Infinity, eta, 2, nu]]]
Sqrt[[Pi]]
1/a14 NIntegrate[
Erfc[u] Exp[(
2 a14 a24 s[2] (-Sqrt[2] Sqrt[a] u + a12 s[2]) -
d (2 a u^2 - 2 Sqrt[2] Sqrt[a] a12 u s[2] + a12^2 s[2]^2) +
a14^2 (2 u^2 - b s[2]^2))/(
2 a14^2) + ((Sqrt[2] Sqrt[
a] (-a14 a34 + a13 d) u + (-a14^2 a23 + a13 a14 a24 +
a12 a14 a34 - a12 a13 d) s[2]) s[3])/
a14^2 - ((-2 a13 a14 a34 + a14^2 c + a13^2 d) s[3]^2)/(
2 a14^2)], u, 0, Infinity, s[2], 0,
Sqrt[2] Sqrt[a]/a12 u, s[3], 0, (Sqrt[2 a] u - a12 s[2])/a13]
Pi/Sqrt[2 dd]
NIntegrate[
Erfc[u] Exp[(A00 u^2 + A10 u s[2] + A11 s[2]^2)/(
2 (dd))] (Erf[(B1 u + B2 s[2])/( a13 Sqrt[2 dd])] +
Erf[(C1 u + C2 s[2])/( a14^1 Sqrt[2 dd])]), u, 0,
Infinity, s[2], 0, Sqrt[2] Sqrt[a]/a12 u]
Now, I will provide the result. Note that the only assumptions on the underlying matrix $bf A$ are that it is symmetric and that its elements are non-negative. Firstly let us define:
begineqnarray
&&mathfrak J^(1,1)(a,b,c)= frac1pi^2cdot left(right.\
&&left.
-frac18
sumlimits_i=1^4 sumlimits_j=1^4
(-1)^j-1+lfloor fraci-12 rfloor
%
mathfrak F^(1,fracsqrt1+2 a^2+b^2 - sqrt2 asqrt1+b^2)_fraci sqrtb^2 c^2+b^2+1 (-1)^leftlfloor
fracj-12rightrfloor +i b c (-1)^jsqrtb^2+1,-fracb
(-1)^i+i (-1)^leftlceil fraci-12rightrceil
sqrtb^2+1
%
right. \
&&left.
right)quad (ii)
endeqnarray
where $mathfrak F^(A,B)_a,b$ is related to di-logarithms and is defined in An integral involving a Gaussian, error functions and the Owen's T function. .
Then we define another function as follows:
beginequation
bar mathfrak J^(1,1)(a,b,c):= fracpi2 arctanleft[ fracsqrt2 a csqrt2 a+b^2(1+c^2)right] - fracpi2 arctanleft[ cright] - 2 pi^2 mathfrak J^(1,1)(frac1sqrt2 a,fracbsqrt2 a,c)
endequation
and then the following quantities that depend on the underlying matrix. We have:
begineqnarray
delta&:=& a_3,3 a_4,1^2 - 2 a_3,1 a_3,4 a_4,1 + a_4,4 a_3,1^2\
W&:=&left(a_3,3 a_4,4-a_3,4^2right) a_1,2^2+2 a_1,4 (a_2,3
a_3,4-a_2,4 a_3,3) a_1,2+2 a_1,3 (a_2,4 a_3,4-a_2,3 a_4,4)
a_1,2+a_1,4^2 left(a_2,2 a_3,3-a_2,3^2right)+2 a_1,3 a_1,4
(a_2,3 a_2,4-a_2,2 a_3,4)+a_1,3^2 left(a_2,2
a_4,4-a_2,4^2right)\
W_1&:=&2 sqrta_1,1 left(a_1,4 (a_2,4 a_3,3-a_2,3 a_3,4)+a_1,3 (a_2,3
a_4,4-a_2,4 a_3,4)+a_1,2 left(a_3,4^2-a_3,3 a_4,4right)right)\
%
v_1&:=&frac1a_4,1 sqrtdelta left( sqrta_1,1(a_3,4 a_4,1 - a_3,1 a_4,4),-a_2,4 a_3,1 a_4,1 + a_2,3 a_4,1^2+a_2,1(-a_3,4 a_4,1+a_3,1a_4,4)right)\
v_2&:=&-frac1a_3,1 sqrtdelta left(sqrta_1,1(a_3,4 a_3,1 - a_4,1 a_3,3),-a_3,1 a_3,2 a_4,1 +a_2,4 a_3,1^2 + a_2,1(-a_3,4 a_3,1+a_4,1a_3,3) right)\
%
left( A, B right)&:=& frac1delta left( W,W_1 right)\
left( bf a_1,bf a_2 right)&:=& frac1sqrtA left(v_1(2),v_2(2) right)\
bf b_1&:=& sqrt2 v_1(1) - fracBsqrt2 A v_1(2)\
bf b_2&:=& sqrt2 v_2(1) - fracBsqrt2 A v_2(2)\
x&:=& fracsqrta_1,1a_2,1
endeqnarray
Then the result reads:
begineqnarray
&&P=frac1det(bf A) left(right.\
%
&& bar mathfrak J^(1,1)left( fracdet(bf A)W,fracBsqrt2 A,bf a_2+fracsqrt2 A bf b_2Bright) -
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracB+2 A xsqrt2 A,bf a_2+fracsqrt2 A bf b_2B+2 A xright)+\
&&!!!!!!!!!!!!!!!!!!!!!!!!!!! bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_2sqrt1+bf a_2^2,bf a_2+fracB(1+bf a_2^2)sqrt2 Abf b_2right) -
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_2sqrt1+bf a_2^2,bf a_2+frac(B+2 A x)(1+bf a_2^2)sqrt2 Abf b_2right)+\
%
&& -bar mathfrak J^(1,1)left( fracdet(bf A)W,fracBsqrt2 A,bf a_1+fracsqrt2 A bf b_1Bright) +
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracB+2 A xsqrt2 A,bf a_1+fracsqrt2 A bf b_1B+2 A xright)+\
&&!!!!!!!!!!!!!!!!!!!!!!!!!!! -bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_1sqrt1+bf a_1^2,bf a_1+fracB(1+bf a_1^2)sqrt2 Abf b_1right) +
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_1sqrt1+bf a_1^2,bf a_1+frac(B+2 A x)(1+bf a_1^2)sqrt2 Abf b_1right)\
%
&&left.right)
endeqnarray
I can provide a code for testing the above expression if anyone is interested.
Now, in the particular case when all the diagonal elements of the matrix $bf A$ are equal unity and all the cross diagonal terms are equal to $rho$ where $0 le rho le 1$ then the result reads:
begineqnarray
&&P=\
&&frac2 pi ^3/2sqrt(1-rho )^3 (3 rho +1)
left( fracpi -3 arctanleft(sqrtfrac3 rho +1rho +1right)2 sqrtpi +6sqrtpi mathfrak J^(1,1)left( fracsqrtfrac32 rho sqrt(1-rho ) (3 rho +1),fracsqrt1-rho sqrt2 sqrt(1-rho ) (3 rho +1),sqrt3right)right)
endeqnarray
Below I plot the quantity $P$ as a function of $rho$. Note that the value $P(rho=0) = pi^2/4 simeq 2.4674$ as it is.
answered Jul 13 '17 at 16:19
PrzemoPrzemo
4,40811031
$endgroup$
add a comment |
$begingroup$
Let us compute the result in case $n=2$. Here the matrix reads $A=left(beginarrayrra & c\c& bendarrayright)$ .Therefore we have:
begineqnarray
P&=& intlimits_mathbb R_+^2 expleft-frac12left[sqrta(s_1+fracca s_2)right]^2 -frac12 fracb a-c^2a s_2^2right ds_1 ds_2\
&=&frac1sqrta sqrtfracpi2 intlimits_0^infty erfcleft(fraccsqrta fracs_2sqrt2 right)expleft-frac12(fracb a-c^2a)s_2^2 rightds_2\
&=&sqrtfracpi2 frac1sqrtb a-c^2 intlimits_0^infty erfc(fraccsqrtb a-c^2 fracs_2sqrt2) e^-frac12 s_2^2 ds_2\
&=& sqrtfracpi2 frac1sqrtb a-c^2 left(
sqrtfracpi2- sqrtfrac2pi arctan(fraccsqrtb a-c^2)right)\
&=& frac1sqrtb a-c^2 arctan(fracsqrtb a-c^2c)
endeqnarray
In the top line we completed the first integration variable to a square and in the second line we integrated over that variable. In the third line we changed variables accordingly . In the fourth line we integrated over the second variable by writing $erfc() = 1- erf()$ and then expanding the error function in a Taylor series and integrating term by term and finally in the last line we simplified the result.
Now, by doing similar calculations we obtained the following result in case $n=3$. Here $A=left(beginarrayrrra & a_12 & a_13\a_12& b&a_23\a_13&a_23&cendarrayright)$.
Firstly we have:
begineqnarray
&&vecs^(T).(A.vecs) = \
&&left(sqrta ( s_1 + fraca_1,2 s_2 + a_1,3 s_3a )right)^2 +
left( b- fraca_1,2^2aright) s_2^2 + left(c-fraca_1,3^2aright) s_3^2 + 2 left(a_2,3-fraca_1,2 a_1,3 aright) s_2 s_3
endeqnarray
Therefore integrating over $s_1$ gives:
begineqnarray
&&P=sqrtfracpi 2 frac1sqrta cdot \
&&intlimits_bf R^2 texterfcleft(fraca_1,2 s_2+a_1,3 s_3sqrt2 sqrtaright) cdot \
&&exp left[
-frac12 left(s_2^2 left(b-fraca_1,2^2aright)+2
s_2 s_3 left(a_2,3-fraca_1,2 a_1,3aright)+s_3^2 left(c-fraca_1,3^2aright)right)
right]
ds_2 ds_3=\
&&
fracsqrtpi a_1,2 intlimits_0^infty
texterfc(u) cdot expleft[-frac12u^2 (frac2 a ba_1,2^2 - 2)right]\
&&
intlimits_0^fracsqrt2 aa_1,3 u
exp left[-frac12 left(s_3 ufrac2 sqrt2 sqrta a_1,2 left(a_2,3-fracb
a_1,3a_1,2right)+
s_3^2fraca_1,3 a_1,2 left(fraca_1,3 ba_1,2+fraca_1,2 ca_1,3-2 a_2,3right)right)right] ds_3 du
endeqnarray
Now it is clear that we can do the integral over $s_3$ in the sense that we can express it through a difference of error functions.Denote $delta:=-2 a_1,2 a_1,3 a_2,3 +a_1,3^2 b +a_1,2^2 c$. Then we have
begineqnarray
&&P=fracpisqrt2sqrtdelta cdotintlimits_0^infty erfc(u) left( erfleft[fracsqrta(-a_1,3 a_2,3+a_1,2 c)a_1,3 sqrtdelta u right] - erfleft[ fracsqrta(a_1,2 a_2,3-a_1,3 b)a_1,2 sqrtdelta u right]right) e^-fracdet(A) delta u^2 du=\
&&fracpisqrt2 det(A)cdot \
&&
intlimits_0^infty erfcleft(u sqrtfracdeltadet(A)right)e^-u^2cdot \
&&left(-erfc(sqrta frac(-a_13a_23+a_12 c)a_13 sqrtdet(A) u)+erfc(sqrta frac(a_12a_23-a_13 b)a_12 sqrtdet(A) u)right) du \
&&=sqrtfracpi2 det(A)\
left[right.\
&&-arctanleft(fraca_13 sqrtdet(A)sqrta(-a_13a_23+a_12 c)right)+
arctanleft(fracsqrtc sqrtdet(A)-a_13 a_23 + a_12 cright)
\
&&+arctanleft(fraca_12 sqrtdet(A)sqrta (a_12 a_23 - a_13 b)right)-arctanleft(fracsqrtb sqrtdet(A)a_12 a_23 - a_13 bright)
left. right]\
&&=sqrtfracpi2 det(A)\
&&left[right.\
&&left.
arctanleft(frac(a_1,3-sqrta_1,1a_3,3)(a_1,3a_2,3-a_1,2a_3,3)sqrta_1,1 (a_1,3a_2,3-a_1,2a_3,3)^2+a_1,3 sqrta_3,3 det(A) sqrtdet(A)right)+right.\
&&left.
arctanleft(frac(a_1,2-sqrta_1,1a_2,2)(a_1,2a_2,3-a_1,3a_2,2)sqrta_1,1 (a_1,2a_2,3-a_1,3a_2,2)^2+a_1,2 sqrta_2,2 det(A) sqrtdet(A)right)
right]
endeqnarray
where in the last line we used An integral involving error functions and a Gaussian .
I also include a Mathematica code snippet that verifies all the steps involved:
(*3d*)
A =.; B =.; CC =.; A12 =.; A23 =.; A13 =.;
For[DDet = 0, True, ,
A, B, CC, A12, A23, A13 =
RandomReal[0, 1, 6, WorkingPrecision -> 50];
DDet = Det[A, A12, A13, A12, B, A23, A13, A23, CC];
If[DDet > 0, Break[]];
];
a = Sqrt[(-2 A12 A13 A23 + A13^2 B + A12^2 CC)/DDet];
b1, b2 = ( Sqrt[A] (-A13 A23 + A12 CC))/ Sqrt[DDet], (
Sqrt[A] (A12 A23 - A13 B))/ Sqrt[DDet];
AA1, AA2 = 2 Sqrt[2] Sqrt[
A] (( A23 A12 - A13 B)/A12^2), (-2 A12 A13 A23 + A13^2 B +
A12^2 CC)/A12^2;
DDet, a, b1, b2;
NIntegrate[
Exp[-1/2 (A s1^2 + B s2^2 + CC s3^2 + 2 A12 s1 s2 + 2 A23 s2 s3 +
2 A13 s1 s3)], s1, 0, Infinity, s2, 0, Infinity, s3, 0,
Infinity]
NIntegrate[
Exp[-1/2 ((Sqrt[A] (s1 + (A12 s2 + A13 s3)/A))^2 + (B -
A12^2/A) s2^2 + (CC - A13^2/A) s3^2 +
2 (A23 - A12 A13/A) s2 s3)], s1, 0, Infinity, s2, 0,
Infinity, s3, 0, Infinity]
NIntegrate[
1/Sqrt[A] Sqrt[
Pi/2] Erfc[(A12 s2 + A13 s3)/
Sqrt[2 A]] Exp[-1/
2 ((B - A12^2/A) s2^2 + (CC - A13^2/A) s3^2 +
2 (A23 - A12 A13/A) s2 s3)], s2, 0, Infinity, s3, 0,
Infinity]
Sqrt[Pi]/A12 NIntegrate[
Erfc[u] Exp[-1/
2 ( A13/A12 (-2 A23 + (A13 B)/A12 + CC A12/A13) s3^2 + (
2 Sqrt[2] Sqrt[A] )/
A12 ( A23 - ( A13 B)/A12) s3 u + (-2 + (2 A B)/
A12^2) u^2)], u, 0, Infinity, s3, 0, Sqrt[2 A]/A13 u]
Sqrt[Pi]/A12 NIntegrate[
Erfc[u] Exp[-1/2 (Sqrt[AA2] s3 + u/2 AA1/Sqrt[AA2])^2] Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity, s3, 0, Sqrt[2 A]/A13 u]
Sqrt[Pi]/(A12 Sqrt[AA2])
NIntegrate[
Erfc[u] Exp[-1/2 (s3)^2] Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity, s3,
u/2 AA1/Sqrt[AA2], ((A13 AA1 + 2 AA2 Sqrt[2] Sqrt[A]) u)/(
2 A13 Sqrt[AA2])]
Sqrt[Pi]/(A12 Sqrt[AA2]) Sqrt[[Pi]/2]
NIntegrate[
Erfc[u] (
Erf[(A13 AA1 + 2 AA2 Sqrt[2] Sqrt[A])/(2 A13 Sqrt[2] Sqrt[AA2])
u] - Erf[AA1/(2 Sqrt[2] Sqrt[AA2]) u]) Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity]
Pi/Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC] Sqrt[1/2]
NIntegrate[
Erfc[u] (
Erf[( Sqrt[A] (-A13 A23 + A12 CC) u)/(
A13 Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC])] -
Erf[(Sqrt[A] (A12 A23 - A13 B) u)/(
A12 Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC])]) Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity]
Pi/ Sqrt[-2 A12 A13 A23 + A13^2 B +
A12^2 CC] Sqrt[1/2] a NIntegrate[
Erfc[a u] (
Erf[( Sqrt[A] (-A13 A23 + A12 CC) u)/(A13 Sqrt[DDet])] -
Erf[(Sqrt[A] (A12 A23 - A13 B) u)/(A12 Sqrt[DDet])]) Exp[-
u^2], u, 0, Infinity]
Pi/Sqrt[2 DDet] NIntegrate[(Erfc[u a]) Exp[-u^2] (Erf[b1/A13 u] -
Erf[b2/A12 u]), u, 0, Infinity]
Sqrt[Pi]/Sqrt[
2 DDet] (ArcTan[ Sqrt[A]/A13 (-A13 A23 + A12 CC)/ Sqrt[DDet]] -
ArcTan[1/ Sqrt[CC] (-A13 A23 + A12 CC)/ Sqrt[DDet]] -
ArcTan[ Sqrt[A]/A12 (A12 A23 - A13 B)/ Sqrt[DDet]] +
ArcTan[ 1/Sqrt[B] (A12 A23 - A13 B)/ Sqrt[DDet]])
-(Sqrt[Pi]/
Sqrt[2 DDet]) (ArcTan[(A13 Sqrt[DDet])/(
Sqrt[A] (-A13 A23 + A12 CC))] -
ArcTan[(Sqrt[CC] Sqrt[DDet])/(-A13 A23 + A12 CC)] -
ArcTan[(A12 Sqrt[DDet])/(Sqrt[A] (A12 A23 - A13 B))] +
ArcTan[(Sqrt[B] Sqrt[DDet])/(A12 A23 - A13 B)])
Sqrt[Pi]/Sqrt[
2 DDet] (ArcTan[((A13 - Sqrt[A] Sqrt[CC]) (A13 A23 - A12 CC) Sqrt[
DDet])/(Sqrt[A] (A13 A23 - A12 CC)^2 + A13 Sqrt[CC] DDet)] +
ArcTan[((A12 - Sqrt[A] Sqrt[B]) (A12 A23 - A13 B) Sqrt[DDet])/(
Sqrt[A] (A12 A23 - A13 B)^2 + A12 Sqrt[B] DDet)])
Update: Now let us take a look at the $n=4$ case. In here:
beginequation
bf A=left(
beginarrayrrrr
a & a_1,2 & a_1,3 & a_1,4 \
a_1,2 & b & a_2,3 & a_2,4 \
a_1,3 & a_2,3 & c & a_3,4 \
a_1,4 & a_2,4 & a_3,4 & d
endarray
right)
endequation
then by doing basically the same computations as a above we managed to reduce the integral in question to a following two dimensional integral. We have:
begineqnarray
&&P= \
&&!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!fracpisqrt2 delta intlimits_0^infty intlimits_0^fracsqrt2 aa_1,2 u
erfc[u] cdot expleft[fracmathfrak A_0,0 u^2 + mathfrak A_1,0 u s_2 +mathfrak A_1,1 s_2^22 deltaright] cdot left( erf[fracmathfrak B_1 u + mathfrak B_2 s_2a_1,3 sqrt2 delta] + erf[fracmathfrak C_1 u + mathfrak C_2 s_2a_1,4 sqrt2 delta]right) d s_2 du=\
&&frac2 imath pi^3/2sqrtmathfrak A_1,1
intlimits_0^infty erfc[u] expfrac4 mathfrak A_0,0 mathfrak A_1,1 - mathfrak A_1,0^28delta mathfrak A_1,1 u^2cdot\
&&
left[right.\
&&left.
left.Tleft(frac(mathfrak A_1,0+xi) u2imath sqrtmathfrak A_1,1 delta,
fracimath mathfrak B_2a_1,3 sqrtmathfrak A_1,1,fracu(2mathfrak A_1,1 mathfrak B_1 - mathfrak A_1,0 mathfrak B_2)2sqrtdelta a_1,3 mathfrak A_1,1right)right|_frac2mathfrak A_1,1 sqrt2 aa_1,2^0 +right.\
&&left.
left.Tleft(frac(mathfrak A_1,0+xi) u2imath sqrtmathfrak A_1,1 delta,
fracimath mathfrak C_2a_1,3 sqrtmathfrak A_1,1,fracu(2mathfrak A_1,1 mathfrak C_1 - mathfrak A_1,0 mathfrak C_2)2sqrtdelta a_1,3 mathfrak A_1,1right)right|_frac2mathfrak A_1,1 sqrt2 aa_1,2^0 +right.\
&&left.
right] du quad (i)
endeqnarray
where $T(cdot,cdot,cdot)$ is the generalized Owen's T function Generalized Owen's T function and
begineqnarray
delta&:=&a_1,3(a_1,3 d-a_1,4 a_3,4) + a_1,4(a_1,4 c- a_1,3 a_3,4)\
mathfrak A_0,0&:=&2 a left(a_3,4^2-c dright)+2 a_1,4 (a_1,4 c-a_1,3 a_3,4)+2 a_1,3 (a_1,3 d-a_1,4 a_3,4)\
mathfrak A_1,0&:=&2 sqrt2 sqrta left(a_1,2 left(c d-a_3,4^2right)+a_1,3
(a_2,4 a_3,4-a_2,3 d)+a_1,4 (a_2,3 a_3,4-a_2,4 c)right)\
mathfrak A_1,1&:=&a_1,2^2 left(a_3,4^2-c dright)+2 a_1,2 a_1,3 (a_2,3 d-a_2,4 a_3,4)+2 a_1,2 a_1,4
(a_2,4 c-a_2,3 a_3,4)+a_1,3^2 left(a_2,4^2-b dright)+2 a_1,3 a_1,4 (a_3,4 b-a_2,3 a_2,4)+a_1,4^2 left(a_2,3^2-b cright)\
hline\
mathfrak B_1&:=&sqrt2 sqrta (a_1,4
c-a_1,3 a_3,4)\
mathfrak B_2&:=&a_1,2 (a_1,3 a_3,4-a_1,4 c)+a_1,3 (a_1,4 a_2,3-a_1,3 a_2,4)\
mathfrak C_1&:=&sqrt2 sqrta (a_1,3 d-a_1,4 a_3,4)\
mathfrak C_2&:=&a_1,2 (a_1,4
a_3,4-a_1,3 d)+a_1,4 (a_1,3 a_2,4-a_1,4 a_2,3)
endeqnarray
nu = 4; Clear[T]; Clear[a]; x =.;
(*a0.dat, a1.dat or a2.dat*)
mat = << "a0.dat";
a, b, c, d, a12, a13, a14, a23, a24, a34 = mat[[1, 1]],
mat[[2, 2]], mat[[3, 3]], mat[[4, 4]], mat[[1, 2]], mat[[1, 3]],
mat[[1, 4]], mat[[2, 3]], mat[[2, 4]], mat[[3, 4]];
dd, A00, A10,
A11 = -2 a13 a14 a34 + a14^2 c + a13^2 d, -4 a13 a14 a34 +
2 a a34^2 + 2 a14^2 c + 2 a13^2 d - 2 a c d,
2 Sqrt[2] Sqrt[a] a14 a23 a34 + 2 Sqrt[2] Sqrt[a] a13 a24 a34 -
2 Sqrt[2] Sqrt[a] a12 a34^2 - 2 Sqrt[2] Sqrt[a] a14 a24 c -
2 Sqrt[2] Sqrt[a] a13 a23 d + 2 Sqrt[2] Sqrt[a] a12 c d,
a14^2 a23^2 - 2 a13 a14 a23 a24 + a13^2 a24^2 -
2 a12 a14 a23 a34 - 2 a12 a13 a24 a34 + a12^2 a34^2 +
2 a13 a14 a34 b + 2 a12 a14 a24 c - a14^2 b c + 2 a12 a13 a23 d -
a13^2 b d - a12^2 c d;
B1, B2, C1,
C2 = Sqrt[2] Sqrt[
a] (-a13 a34 + a14 c), (a13 a14 a23 - a13^2 a24 + a12 a13 a34 -
a12 a14 c),
Sqrt[2] Sqrt[
a] (-a14 a34 + a13 d), (-a14^2 a23 + a13 a14 a24 + a12 a14 a34 -
a12 a13 d);
NIntegrate[
Exp[-1/2 Sum[mat[[i, j]] s[i] s[j], i, 1, nu, j, 1, nu]],
Evaluate[Sequence @@ Table[s[eta], 0, Infinity, eta, 1, nu]]]
Sqrt[[Pi]/(2 a)]
NIntegrate[
Erfc[(a12 s[2] + a13 s[3] + a14 s[4])/Sqrt[
2 a]] Exp[-1/
2 ((-(a12^2/a) + b) s[2]^2 + (-(a13^2/a) + c) s[
3]^2 + (-(a14^2/a) + d) s[4]^2 +
2 (-(( a13 a14)/a) + a34) s[3] s[4] +
2 (-(( a12 a13)/a) + a23) s[2] s[3] +
2 (-(( a12 a14)/a) + a24) s[2] s[4])],
Evaluate[Sequence @@ Table[s[eta], 0, Infinity, eta, 2, nu]]]
Sqrt[[Pi]]
1/a14 NIntegrate[
Erfc[u] Exp[(
2 a14 a24 s[2] (-Sqrt[2] Sqrt[a] u + a12 s[2]) -
d (2 a u^2 - 2 Sqrt[2] Sqrt[a] a12 u s[2] + a12^2 s[2]^2) +
a14^2 (2 u^2 - b s[2]^2))/(
2 a14^2) + ((Sqrt[2] Sqrt[
a] (-a14 a34 + a13 d) u + (-a14^2 a23 + a13 a14 a24 +
a12 a14 a34 - a12 a13 d) s[2]) s[3])/
a14^2 - ((-2 a13 a14 a34 + a14^2 c + a13^2 d) s[3]^2)/(
2 a14^2)], u, 0, Infinity, s[2], 0,
Sqrt[2] Sqrt[a]/a12 u, s[3], 0, (Sqrt[2 a] u - a12 s[2])/a13]
Pi/Sqrt[2 dd]
NIntegrate[
Erfc[u] Exp[(A00 u^2 + A10 u s[2] + A11 s[2]^2)/(
2 (dd))] (Erf[(B1 u + B2 s[2])/( a13 Sqrt[2 dd])] +
Erf[(C1 u + C2 s[2])/( a14^1 Sqrt[2 dd])]), u, 0,
Infinity, s[2], 0, Sqrt[2] Sqrt[a]/a12 u]
Now, I will provide the result. Note that the only assumptions on the underlying matrix $bf A$ are that it is symmetric and that its elements are non-negative. Firstly let us define:
begineqnarray
&&mathfrak J^(1,1)(a,b,c)= frac1pi^2cdot left(right.\
&&left.
-frac18
sumlimits_i=1^4 sumlimits_j=1^4
(-1)^j-1+lfloor fraci-12 rfloor
%
mathfrak F^(1,fracsqrt1+2 a^2+b^2 - sqrt2 asqrt1+b^2)_fraci sqrtb^2 c^2+b^2+1 (-1)^leftlfloor
fracj-12rightrfloor +i b c (-1)^jsqrtb^2+1,-fracb
(-1)^i+i (-1)^leftlceil fraci-12rightrceil
sqrtb^2+1
%
right. \
&&left.
right)quad (ii)
endeqnarray
where $mathfrak F^(A,B)_a,b$ is related to di-logarithms and is defined in An integral involving a Gaussian, error functions and the Owen's T function. .
Then we define another function as follows:
beginequation
bar mathfrak J^(1,1)(a,b,c):= fracpi2 arctanleft[ fracsqrt2 a csqrt2 a+b^2(1+c^2)right] - fracpi2 arctanleft[ cright] - 2 pi^2 mathfrak J^(1,1)(frac1sqrt2 a,fracbsqrt2 a,c)
endequation
and then the following quantities that depend on the underlying matrix. We have:
begineqnarray
delta&:=& a_3,3 a_4,1^2 - 2 a_3,1 a_3,4 a_4,1 + a_4,4 a_3,1^2\
W&:=&left(a_3,3 a_4,4-a_3,4^2right) a_1,2^2+2 a_1,4 (a_2,3
a_3,4-a_2,4 a_3,3) a_1,2+2 a_1,3 (a_2,4 a_3,4-a_2,3 a_4,4)
a_1,2+a_1,4^2 left(a_2,2 a_3,3-a_2,3^2right)+2 a_1,3 a_1,4
(a_2,3 a_2,4-a_2,2 a_3,4)+a_1,3^2 left(a_2,2
a_4,4-a_2,4^2right)\
W_1&:=&2 sqrta_1,1 left(a_1,4 (a_2,4 a_3,3-a_2,3 a_3,4)+a_1,3 (a_2,3
a_4,4-a_2,4 a_3,4)+a_1,2 left(a_3,4^2-a_3,3 a_4,4right)right)\
%
v_1&:=&frac1a_4,1 sqrtdelta left( sqrta_1,1(a_3,4 a_4,1 - a_3,1 a_4,4),-a_2,4 a_3,1 a_4,1 + a_2,3 a_4,1^2+a_2,1(-a_3,4 a_4,1+a_3,1a_4,4)right)\
v_2&:=&-frac1a_3,1 sqrtdelta left(sqrta_1,1(a_3,4 a_3,1 - a_4,1 a_3,3),-a_3,1 a_3,2 a_4,1 +a_2,4 a_3,1^2 + a_2,1(-a_3,4 a_3,1+a_4,1a_3,3) right)\
%
left( A, B right)&:=& frac1delta left( W,W_1 right)\
left( bf a_1,bf a_2 right)&:=& frac1sqrtA left(v_1(2),v_2(2) right)\
bf b_1&:=& sqrt2 v_1(1) - fracBsqrt2 A v_1(2)\
bf b_2&:=& sqrt2 v_2(1) - fracBsqrt2 A v_2(2)\
x&:=& fracsqrta_1,1a_2,1
endeqnarray
Then the result reads:
begineqnarray
&&P=frac1det(bf A) left(right.\
%
&& bar mathfrak J^(1,1)left( fracdet(bf A)W,fracBsqrt2 A,bf a_2+fracsqrt2 A bf b_2Bright) -
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracB+2 A xsqrt2 A,bf a_2+fracsqrt2 A bf b_2B+2 A xright)+\
&&!!!!!!!!!!!!!!!!!!!!!!!!!!! bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_2sqrt1+bf a_2^2,bf a_2+fracB(1+bf a_2^2)sqrt2 Abf b_2right) -
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_2sqrt1+bf a_2^2,bf a_2+frac(B+2 A x)(1+bf a_2^2)sqrt2 Abf b_2right)+\
%
&& -bar mathfrak J^(1,1)left( fracdet(bf A)W,fracBsqrt2 A,bf a_1+fracsqrt2 A bf b_1Bright) +
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracB+2 A xsqrt2 A,bf a_1+fracsqrt2 A bf b_1B+2 A xright)+\
&&!!!!!!!!!!!!!!!!!!!!!!!!!!! -bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_1sqrt1+bf a_1^2,bf a_1+fracB(1+bf a_1^2)sqrt2 Abf b_1right) +
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_1sqrt1+bf a_1^2,bf a_1+frac(B+2 A x)(1+bf a_1^2)sqrt2 Abf b_1right)\
%
&&left.right)
endeqnarray
I can provide a code for testing the above expression if anyone is interested.
Now, in the particular case when all the diagonal elements of the matrix $bf A$ are equal unity and all the cross diagonal terms are equal to $rho$ where $0 le rho le 1$ then the result reads:
begineqnarray
&&P=\
&&frac2 pi ^3/2sqrt(1-rho )^3 (3 rho +1)
left( fracpi -3 arctanleft(sqrtfrac3 rho +1rho +1right)2 sqrtpi +6sqrtpi mathfrak J^(1,1)left( fracsqrtfrac32 rho sqrt(1-rho ) (3 rho +1),fracsqrt1-rho sqrt2 sqrt(1-rho ) (3 rho +1),sqrt3right)right)
endeqnarray
Below I plot the quantity $P$ as a function of $rho$. Note that the value $P(rho=0) = pi^2/4 simeq 2.4674$ as it is.
answered Jul 13 '17 at 16:19
PrzemoPrzemo
4,40811031
$endgroup$
Let us compute the result in case $n=2$. Here the matrix reads $A=left(beginarrayrra & c\c& bendarrayright)$ .Therefore we have:
begineqnarray
P&=& intlimits_mathbb R_+^2 expleft-frac12left[sqrta(s_1+fracca s_2)right]^2 -frac12 fracb a-c^2a s_2^2right ds_1 ds_2\
&=&frac1sqrta sqrtfracpi2 intlimits_0^infty erfcleft(fraccsqrta fracs_2sqrt2 right)expleft-frac12(fracb a-c^2a)s_2^2 rightds_2\
&=&sqrtfracpi2 frac1sqrtb a-c^2 intlimits_0^infty erfc(fraccsqrtb a-c^2 fracs_2sqrt2) e^-frac12 s_2^2 ds_2\
&=& sqrtfracpi2 frac1sqrtb a-c^2 left(
sqrtfracpi2- sqrtfrac2pi arctan(fraccsqrtb a-c^2)right)\
&=& frac1sqrtb a-c^2 arctan(fracsqrtb a-c^2c)
endeqnarray
In the top line we completed the first integration variable to a square and in the second line we integrated over that variable. In the third line we changed variables accordingly . In the fourth line we integrated over the second variable by writing $erfc() = 1- erf()$ and then expanding the error function in a Taylor series and integrating term by term and finally in the last line we simplified the result.
Now, by doing similar calculations we obtained the following result in case $n=3$. Here $A=left(beginarrayrrra & a_12 & a_13\a_12& b&a_23\a_13&a_23&cendarrayright)$.
Firstly we have:
begineqnarray
&&vecs^(T).(A.vecs) = \
&&left(sqrta ( s_1 + fraca_1,2 s_2 + a_1,3 s_3a )right)^2 +
left( b- fraca_1,2^2aright) s_2^2 + left(c-fraca_1,3^2aright) s_3^2 + 2 left(a_2,3-fraca_1,2 a_1,3 aright) s_2 s_3
endeqnarray
Therefore integrating over $s_1$ gives:
begineqnarray
&&P=sqrtfracpi 2 frac1sqrta cdot \
&&intlimits_bf R^2 texterfcleft(fraca_1,2 s_2+a_1,3 s_3sqrt2 sqrtaright) cdot \
&&exp left[
-frac12 left(s_2^2 left(b-fraca_1,2^2aright)+2
s_2 s_3 left(a_2,3-fraca_1,2 a_1,3aright)+s_3^2 left(c-fraca_1,3^2aright)right)
right]
ds_2 ds_3=\
&&
fracsqrtpi a_1,2 intlimits_0^infty
texterfc(u) cdot expleft[-frac12u^2 (frac2 a ba_1,2^2 - 2)right]\
&&
intlimits_0^fracsqrt2 aa_1,3 u
exp left[-frac12 left(s_3 ufrac2 sqrt2 sqrta a_1,2 left(a_2,3-fracb
a_1,3a_1,2right)+
s_3^2fraca_1,3 a_1,2 left(fraca_1,3 ba_1,2+fraca_1,2 ca_1,3-2 a_2,3right)right)right] ds_3 du
endeqnarray
Now it is clear that we can do the integral over $s_3$ in the sense that we can express it through a difference of error functions.Denote $delta:=-2 a_1,2 a_1,3 a_2,3 +a_1,3^2 b +a_1,2^2 c$. Then we have
begineqnarray
&&P=fracpisqrt2sqrtdelta cdotintlimits_0^infty erfc(u) left( erfleft[fracsqrta(-a_1,3 a_2,3+a_1,2 c)a_1,3 sqrtdelta u right] - erfleft[ fracsqrta(a_1,2 a_2,3-a_1,3 b)a_1,2 sqrtdelta u right]right) e^-fracdet(A) delta u^2 du=\
&&fracpisqrt2 det(A)cdot \
&&
intlimits_0^infty erfcleft(u sqrtfracdeltadet(A)right)e^-u^2cdot \
&&left(-erfc(sqrta frac(-a_13a_23+a_12 c)a_13 sqrtdet(A) u)+erfc(sqrta frac(a_12a_23-a_13 b)a_12 sqrtdet(A) u)right) du \
&&=sqrtfracpi2 det(A)\
left[right.\
&&-arctanleft(fraca_13 sqrtdet(A)sqrta(-a_13a_23+a_12 c)right)+
arctanleft(fracsqrtc sqrtdet(A)-a_13 a_23 + a_12 cright)
\
&&+arctanleft(fraca_12 sqrtdet(A)sqrta (a_12 a_23 - a_13 b)right)-arctanleft(fracsqrtb sqrtdet(A)a_12 a_23 - a_13 bright)
left. right]\
&&=sqrtfracpi2 det(A)\
&&left[right.\
&&left.
arctanleft(frac(a_1,3-sqrta_1,1a_3,3)(a_1,3a_2,3-a_1,2a_3,3)sqrta_1,1 (a_1,3a_2,3-a_1,2a_3,3)^2+a_1,3 sqrta_3,3 det(A) sqrtdet(A)right)+right.\
&&left.
arctanleft(frac(a_1,2-sqrta_1,1a_2,2)(a_1,2a_2,3-a_1,3a_2,2)sqrta_1,1 (a_1,2a_2,3-a_1,3a_2,2)^2+a_1,2 sqrta_2,2 det(A) sqrtdet(A)right)
right]
endeqnarray
where in the last line we used An integral involving error functions and a Gaussian .
I also include a Mathematica code snippet that verifies all the steps involved:
(*3d*)
A =.; B =.; CC =.; A12 =.; A23 =.; A13 =.;
For[DDet = 0, True, ,
A, B, CC, A12, A23, A13 =
RandomReal[0, 1, 6, WorkingPrecision -> 50];
DDet = Det[A, A12, A13, A12, B, A23, A13, A23, CC];
If[DDet > 0, Break[]];
];
a = Sqrt[(-2 A12 A13 A23 + A13^2 B + A12^2 CC)/DDet];
b1, b2 = ( Sqrt[A] (-A13 A23 + A12 CC))/ Sqrt[DDet], (
Sqrt[A] (A12 A23 - A13 B))/ Sqrt[DDet];
AA1, AA2 = 2 Sqrt[2] Sqrt[
A] (( A23 A12 - A13 B)/A12^2), (-2 A12 A13 A23 + A13^2 B +
A12^2 CC)/A12^2;
DDet, a, b1, b2;
NIntegrate[
Exp[-1/2 (A s1^2 + B s2^2 + CC s3^2 + 2 A12 s1 s2 + 2 A23 s2 s3 +
2 A13 s1 s3)], s1, 0, Infinity, s2, 0, Infinity, s3, 0,
Infinity]
NIntegrate[
Exp[-1/2 ((Sqrt[A] (s1 + (A12 s2 + A13 s3)/A))^2 + (B -
A12^2/A) s2^2 + (CC - A13^2/A) s3^2 +
2 (A23 - A12 A13/A) s2 s3)], s1, 0, Infinity, s2, 0,
Infinity, s3, 0, Infinity]
NIntegrate[
1/Sqrt[A] Sqrt[
Pi/2] Erfc[(A12 s2 + A13 s3)/
Sqrt[2 A]] Exp[-1/
2 ((B - A12^2/A) s2^2 + (CC - A13^2/A) s3^2 +
2 (A23 - A12 A13/A) s2 s3)], s2, 0, Infinity, s3, 0,
Infinity]
Sqrt[Pi]/A12 NIntegrate[
Erfc[u] Exp[-1/
2 ( A13/A12 (-2 A23 + (A13 B)/A12 + CC A12/A13) s3^2 + (
2 Sqrt[2] Sqrt[A] )/
A12 ( A23 - ( A13 B)/A12) s3 u + (-2 + (2 A B)/
A12^2) u^2)], u, 0, Infinity, s3, 0, Sqrt[2 A]/A13 u]
Sqrt[Pi]/A12 NIntegrate[
Erfc[u] Exp[-1/2 (Sqrt[AA2] s3 + u/2 AA1/Sqrt[AA2])^2] Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity, s3, 0, Sqrt[2 A]/A13 u]
Sqrt[Pi]/(A12 Sqrt[AA2])
NIntegrate[
Erfc[u] Exp[-1/2 (s3)^2] Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity, s3,
u/2 AA1/Sqrt[AA2], ((A13 AA1 + 2 AA2 Sqrt[2] Sqrt[A]) u)/(
2 A13 Sqrt[AA2])]
Sqrt[Pi]/(A12 Sqrt[AA2]) Sqrt[[Pi]/2]
NIntegrate[
Erfc[u] (
Erf[(A13 AA1 + 2 AA2 Sqrt[2] Sqrt[A])/(2 A13 Sqrt[2] Sqrt[AA2])
u] - Erf[AA1/(2 Sqrt[2] Sqrt[AA2]) u]) Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity]
Pi/Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC] Sqrt[1/2]
NIntegrate[
Erfc[u] (
Erf[( Sqrt[A] (-A13 A23 + A12 CC) u)/(
A13 Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC])] -
Erf[(Sqrt[A] (A12 A23 - A13 B) u)/(
A12 Sqrt[-2 A12 A13 A23 + A13^2 B + A12^2 CC])]) Exp[-((
DDet u^2)/(-2 A12 A13 A23 + A13^2 B + A12^2 CC))], u, 0,
Infinity]
Pi/ Sqrt[-2 A12 A13 A23 + A13^2 B +
A12^2 CC] Sqrt[1/2] a NIntegrate[
Erfc[a u] (
Erf[( Sqrt[A] (-A13 A23 + A12 CC) u)/(A13 Sqrt[DDet])] -
Erf[(Sqrt[A] (A12 A23 - A13 B) u)/(A12 Sqrt[DDet])]) Exp[-
u^2], u, 0, Infinity]
Pi/Sqrt[2 DDet] NIntegrate[(Erfc[u a]) Exp[-u^2] (Erf[b1/A13 u] -
Erf[b2/A12 u]), u, 0, Infinity]
Sqrt[Pi]/Sqrt[
2 DDet] (ArcTan[ Sqrt[A]/A13 (-A13 A23 + A12 CC)/ Sqrt[DDet]] -
ArcTan[1/ Sqrt[CC] (-A13 A23 + A12 CC)/ Sqrt[DDet]] -
ArcTan[ Sqrt[A]/A12 (A12 A23 - A13 B)/ Sqrt[DDet]] +
ArcTan[ 1/Sqrt[B] (A12 A23 - A13 B)/ Sqrt[DDet]])
-(Sqrt[Pi]/
Sqrt[2 DDet]) (ArcTan[(A13 Sqrt[DDet])/(
Sqrt[A] (-A13 A23 + A12 CC))] -
ArcTan[(Sqrt[CC] Sqrt[DDet])/(-A13 A23 + A12 CC)] -
ArcTan[(A12 Sqrt[DDet])/(Sqrt[A] (A12 A23 - A13 B))] +
ArcTan[(Sqrt[B] Sqrt[DDet])/(A12 A23 - A13 B)])
Sqrt[Pi]/Sqrt[
2 DDet] (ArcTan[((A13 - Sqrt[A] Sqrt[CC]) (A13 A23 - A12 CC) Sqrt[
DDet])/(Sqrt[A] (A13 A23 - A12 CC)^2 + A13 Sqrt[CC] DDet)] +
ArcTan[((A12 - Sqrt[A] Sqrt[B]) (A12 A23 - A13 B) Sqrt[DDet])/(
Sqrt[A] (A12 A23 - A13 B)^2 + A12 Sqrt[B] DDet)])
Update: Now let us take a look at the $n=4$ case. In here:
beginequation
bf A=left(
beginarrayrrrr
a & a_1,2 & a_1,3 & a_1,4 \
a_1,2 & b & a_2,3 & a_2,4 \
a_1,3 & a_2,3 & c & a_3,4 \
a_1,4 & a_2,4 & a_3,4 & d
endarray
right)
endequation
then by doing basically the same computations as a above we managed to reduce the integral in question to a following two dimensional integral. We have:
begineqnarray
&&P= \
&&!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!fracpisqrt2 delta intlimits_0^infty intlimits_0^fracsqrt2 aa_1,2 u
erfc[u] cdot expleft[fracmathfrak A_0,0 u^2 + mathfrak A_1,0 u s_2 +mathfrak A_1,1 s_2^22 deltaright] cdot left( erf[fracmathfrak B_1 u + mathfrak B_2 s_2a_1,3 sqrt2 delta] + erf[fracmathfrak C_1 u + mathfrak C_2 s_2a_1,4 sqrt2 delta]right) d s_2 du=\
&&frac2 imath pi^3/2sqrtmathfrak A_1,1
intlimits_0^infty erfc[u] expfrac4 mathfrak A_0,0 mathfrak A_1,1 - mathfrak A_1,0^28delta mathfrak A_1,1 u^2cdot\
&&
left[right.\
&&left.
left.Tleft(frac(mathfrak A_1,0+xi) u2imath sqrtmathfrak A_1,1 delta,
fracimath mathfrak B_2a_1,3 sqrtmathfrak A_1,1,fracu(2mathfrak A_1,1 mathfrak B_1 - mathfrak A_1,0 mathfrak B_2)2sqrtdelta a_1,3 mathfrak A_1,1right)right|_frac2mathfrak A_1,1 sqrt2 aa_1,2^0 +right.\
&&left.
left.Tleft(frac(mathfrak A_1,0+xi) u2imath sqrtmathfrak A_1,1 delta,
fracimath mathfrak C_2a_1,3 sqrtmathfrak A_1,1,fracu(2mathfrak A_1,1 mathfrak C_1 - mathfrak A_1,0 mathfrak C_2)2sqrtdelta a_1,3 mathfrak A_1,1right)right|_frac2mathfrak A_1,1 sqrt2 aa_1,2^0 +right.\
&&left.
right] du quad (i)
endeqnarray
where $T(cdot,cdot,cdot)$ is the generalized Owen's T function Generalized Owen's T function and
begineqnarray
delta&:=&a_1,3(a_1,3 d-a_1,4 a_3,4) + a_1,4(a_1,4 c- a_1,3 a_3,4)\
mathfrak A_0,0&:=&2 a left(a_3,4^2-c dright)+2 a_1,4 (a_1,4 c-a_1,3 a_3,4)+2 a_1,3 (a_1,3 d-a_1,4 a_3,4)\
mathfrak A_1,0&:=&2 sqrt2 sqrta left(a_1,2 left(c d-a_3,4^2right)+a_1,3
(a_2,4 a_3,4-a_2,3 d)+a_1,4 (a_2,3 a_3,4-a_2,4 c)right)\
mathfrak A_1,1&:=&a_1,2^2 left(a_3,4^2-c dright)+2 a_1,2 a_1,3 (a_2,3 d-a_2,4 a_3,4)+2 a_1,2 a_1,4
(a_2,4 c-a_2,3 a_3,4)+a_1,3^2 left(a_2,4^2-b dright)+2 a_1,3 a_1,4 (a_3,4 b-a_2,3 a_2,4)+a_1,4^2 left(a_2,3^2-b cright)\
hline\
mathfrak B_1&:=&sqrt2 sqrta (a_1,4
c-a_1,3 a_3,4)\
mathfrak B_2&:=&a_1,2 (a_1,3 a_3,4-a_1,4 c)+a_1,3 (a_1,4 a_2,3-a_1,3 a_2,4)\
mathfrak C_1&:=&sqrt2 sqrta (a_1,3 d-a_1,4 a_3,4)\
mathfrak C_2&:=&a_1,2 (a_1,4
a_3,4-a_1,3 d)+a_1,4 (a_1,3 a_2,4-a_1,4 a_2,3)
endeqnarray
nu = 4; Clear[T]; Clear[a]; x =.;
(*a0.dat, a1.dat or a2.dat*)
mat = << "a0.dat";
a, b, c, d, a12, a13, a14, a23, a24, a34 = mat[[1, 1]],
mat[[2, 2]], mat[[3, 3]], mat[[4, 4]], mat[[1, 2]], mat[[1, 3]],
mat[[1, 4]], mat[[2, 3]], mat[[2, 4]], mat[[3, 4]];
dd, A00, A10,
A11 = -2 a13 a14 a34 + a14^2 c + a13^2 d, -4 a13 a14 a34 +
2 a a34^2 + 2 a14^2 c + 2 a13^2 d - 2 a c d,
2 Sqrt[2] Sqrt[a] a14 a23 a34 + 2 Sqrt[2] Sqrt[a] a13 a24 a34 -
2 Sqrt[2] Sqrt[a] a12 a34^2 - 2 Sqrt[2] Sqrt[a] a14 a24 c -
2 Sqrt[2] Sqrt[a] a13 a23 d + 2 Sqrt[2] Sqrt[a] a12 c d,
a14^2 a23^2 - 2 a13 a14 a23 a24 + a13^2 a24^2 -
2 a12 a14 a23 a34 - 2 a12 a13 a24 a34 + a12^2 a34^2 +
2 a13 a14 a34 b + 2 a12 a14 a24 c - a14^2 b c + 2 a12 a13 a23 d -
a13^2 b d - a12^2 c d;
B1, B2, C1,
C2 = Sqrt[2] Sqrt[
a] (-a13 a34 + a14 c), (a13 a14 a23 - a13^2 a24 + a12 a13 a34 -
a12 a14 c),
Sqrt[2] Sqrt[
a] (-a14 a34 + a13 d), (-a14^2 a23 + a13 a14 a24 + a12 a14 a34 -
a12 a13 d);
NIntegrate[
Exp[-1/2 Sum[mat[[i, j]] s[i] s[j], i, 1, nu, j, 1, nu]],
Evaluate[Sequence @@ Table[s[eta], 0, Infinity, eta, 1, nu]]]
Sqrt[[Pi]/(2 a)]
NIntegrate[
Erfc[(a12 s[2] + a13 s[3] + a14 s[4])/Sqrt[
2 a]] Exp[-1/
2 ((-(a12^2/a) + b) s[2]^2 + (-(a13^2/a) + c) s[
3]^2 + (-(a14^2/a) + d) s[4]^2 +
2 (-(( a13 a14)/a) + a34) s[3] s[4] +
2 (-(( a12 a13)/a) + a23) s[2] s[3] +
2 (-(( a12 a14)/a) + a24) s[2] s[4])],
Evaluate[Sequence @@ Table[s[eta], 0, Infinity, eta, 2, nu]]]
Sqrt[[Pi]]
1/a14 NIntegrate[
Erfc[u] Exp[(
2 a14 a24 s[2] (-Sqrt[2] Sqrt[a] u + a12 s[2]) -
d (2 a u^2 - 2 Sqrt[2] Sqrt[a] a12 u s[2] + a12^2 s[2]^2) +
a14^2 (2 u^2 - b s[2]^2))/(
2 a14^2) + ((Sqrt[2] Sqrt[
a] (-a14 a34 + a13 d) u + (-a14^2 a23 + a13 a14 a24 +
a12 a14 a34 - a12 a13 d) s[2]) s[3])/
a14^2 - ((-2 a13 a14 a34 + a14^2 c + a13^2 d) s[3]^2)/(
2 a14^2)], u, 0, Infinity, s[2], 0,
Sqrt[2] Sqrt[a]/a12 u, s[3], 0, (Sqrt[2 a] u - a12 s[2])/a13]
Pi/Sqrt[2 dd]
NIntegrate[
Erfc[u] Exp[(A00 u^2 + A10 u s[2] + A11 s[2]^2)/(
2 (dd))] (Erf[(B1 u + B2 s[2])/( a13 Sqrt[2 dd])] +
Erf[(C1 u + C2 s[2])/( a14^1 Sqrt[2 dd])]), u, 0,
Infinity, s[2], 0, Sqrt[2] Sqrt[a]/a12 u]
Now, I will provide the result. Note that the only assumptions on the underlying matrix $bf A$ are that it is symmetric and that its elements are non-negative. Firstly let us define:
begineqnarray
&&mathfrak J^(1,1)(a,b,c)= frac1pi^2cdot left(right.\
&&left.
-frac18
sumlimits_i=1^4 sumlimits_j=1^4
(-1)^j-1+lfloor fraci-12 rfloor
%
mathfrak F^(1,fracsqrt1+2 a^2+b^2 - sqrt2 asqrt1+b^2)_fraci sqrtb^2 c^2+b^2+1 (-1)^leftlfloor
fracj-12rightrfloor +i b c (-1)^jsqrtb^2+1,-fracb
(-1)^i+i (-1)^leftlceil fraci-12rightrceil
sqrtb^2+1
%
right. \
&&left.
right)quad (ii)
endeqnarray
where $mathfrak F^(A,B)_a,b$ is related to di-logarithms and is defined in An integral involving a Gaussian, error functions and the Owen's T function. .
Then we define another function as follows:
beginequation
bar mathfrak J^(1,1)(a,b,c):= fracpi2 arctanleft[ fracsqrt2 a csqrt2 a+b^2(1+c^2)right] - fracpi2 arctanleft[ cright] - 2 pi^2 mathfrak J^(1,1)(frac1sqrt2 a,fracbsqrt2 a,c)
endequation
and then the following quantities that depend on the underlying matrix. We have:
begineqnarray
delta&:=& a_3,3 a_4,1^2 - 2 a_3,1 a_3,4 a_4,1 + a_4,4 a_3,1^2\
W&:=&left(a_3,3 a_4,4-a_3,4^2right) a_1,2^2+2 a_1,4 (a_2,3
a_3,4-a_2,4 a_3,3) a_1,2+2 a_1,3 (a_2,4 a_3,4-a_2,3 a_4,4)
a_1,2+a_1,4^2 left(a_2,2 a_3,3-a_2,3^2right)+2 a_1,3 a_1,4
(a_2,3 a_2,4-a_2,2 a_3,4)+a_1,3^2 left(a_2,2
a_4,4-a_2,4^2right)\
W_1&:=&2 sqrta_1,1 left(a_1,4 (a_2,4 a_3,3-a_2,3 a_3,4)+a_1,3 (a_2,3
a_4,4-a_2,4 a_3,4)+a_1,2 left(a_3,4^2-a_3,3 a_4,4right)right)\
%
v_1&:=&frac1a_4,1 sqrtdelta left( sqrta_1,1(a_3,4 a_4,1 - a_3,1 a_4,4),-a_2,4 a_3,1 a_4,1 + a_2,3 a_4,1^2+a_2,1(-a_3,4 a_4,1+a_3,1a_4,4)right)\
v_2&:=&-frac1a_3,1 sqrtdelta left(sqrta_1,1(a_3,4 a_3,1 - a_4,1 a_3,3),-a_3,1 a_3,2 a_4,1 +a_2,4 a_3,1^2 + a_2,1(-a_3,4 a_3,1+a_4,1a_3,3) right)\
%
left( A, B right)&:=& frac1delta left( W,W_1 right)\
left( bf a_1,bf a_2 right)&:=& frac1sqrtA left(v_1(2),v_2(2) right)\
bf b_1&:=& sqrt2 v_1(1) - fracBsqrt2 A v_1(2)\
bf b_2&:=& sqrt2 v_2(1) - fracBsqrt2 A v_2(2)\
x&:=& fracsqrta_1,1a_2,1
endeqnarray
Then the result reads:
begineqnarray
&&P=frac1det(bf A) left(right.\
%
&& bar mathfrak J^(1,1)left( fracdet(bf A)W,fracBsqrt2 A,bf a_2+fracsqrt2 A bf b_2Bright) -
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracB+2 A xsqrt2 A,bf a_2+fracsqrt2 A bf b_2B+2 A xright)+\
&&!!!!!!!!!!!!!!!!!!!!!!!!!!! bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_2sqrt1+bf a_2^2,bf a_2+fracB(1+bf a_2^2)sqrt2 Abf b_2right) -
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_2sqrt1+bf a_2^2,bf a_2+frac(B+2 A x)(1+bf a_2^2)sqrt2 Abf b_2right)+\
%
&& -bar mathfrak J^(1,1)left( fracdet(bf A)W,fracBsqrt2 A,bf a_1+fracsqrt2 A bf b_1Bright) +
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracB+2 A xsqrt2 A,bf a_1+fracsqrt2 A bf b_1B+2 A xright)+\
&&!!!!!!!!!!!!!!!!!!!!!!!!!!! -bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_1sqrt1+bf a_1^2,bf a_1+fracB(1+bf a_1^2)sqrt2 Abf b_1right) +
bar mathfrak J^(1,1)left( fracdet(bf A)W,fracbf b_1sqrt1+bf a_1^2,bf a_1+frac(B+2 A x)(1+bf a_1^2)sqrt2 Abf b_1right)\
%
&&left.right)
endeqnarray
I can provide a code for testing the above expression if anyone is interested.
Now, in the particular case when all the diagonal elements of the matrix $bf A$ are equal unity and all the cross diagonal terms are equal to $rho$ where $0 le rho le 1$ then the result reads:
begineqnarray
&&P=\
&&frac2 pi ^3/2sqrt(1-rho )^3 (3 rho +1)
left( fracpi -3 arctanleft(sqrtfrac3 rho +1rho +1right)2 sqrtpi +6sqrtpi mathfrak J^(1,1)left( fracsqrtfrac32 rho sqrt(1-rho ) (3 rho +1),fracsqrt1-rho sqrt2 sqrt(1-rho ) (3 rho +1),sqrt3right)right)
endeqnarray
Below I plot the quantity $P$ as a function of $rho$. Note that the value $P(rho=0) = pi^2/4 simeq 2.4674$ as it is.
answered Jul 13 '17 at 16:19
PrzemoPrzemo
4,40811031
edited Mar 14 at 16:42
answered Jul 13 '17 at 16:19
PrzemoPrzemo
4,40811031
answered Jul 13 '17 at 16:19
PrzemoPrzemo
4,40811031
answered Jul 13 '17 at 16:19
PrzemoPrzemo
4,40811031
4,40811031
add a comment |
add a comment |
$begingroup$
Thank you Przemo for your solution to the problem for $n=2, 3$. While I had no trouble following your derivation in 2D, I'm stuck with the derivation of your intermediate step for $n=3$. I mainly tried two approaches:
Completing the square in one variable, say $x$, leaves me with
$$int_mathbbR_+^2 mathrmdymathrmdz expleft(-frac12 fracmathrmdet,A_3mathrmdet,A_2z^2right) expleft(-frac12 fracmathrmdet, A_2a(y-m z)^2right) left[1 - mathrmerfleft(fraca_12y+a_13zsqrt2aright) right] $$
where $A_2=beginpmatrix a & a_12\ & bendpmatrix$, $A_3$ as you defined it, and $m$ is a function of the coefficients of the matrices. However, I do not know how to proceed from there: expanding the error function to do the integral in y, say, is a nightmare due to the constant term in z; I also did not find a way to do a coordinate transform à la $s=a_12y+a_13z$ or something similar.Indeed, your intermediate solution looks more like you were able to complete the square in two of the variables independently; but what happened to the cross-term? I cannot find a factorisation of the exponent that would allow me to complete two integrals over the half-line with only one variable left in the error function yielded by the integral.
Any help / hint would be greatly appreciated! Thank you in advance.
answered Dec 3 '18 at 18:09
workandheatworkandheat
587
$endgroup$
1
$begingroup$
@ workandheat, the key word here is orthant probability. If you search for this in connection with normal distribution, you get many results. But as I see even for cases with four variables people write long papers. So for more variables I doubt that there are exact results.
$endgroup$
– Karl
Dec 3 '18 at 18:23
$begingroup$
Thanks @Karl for the hint to "orthant" probability; coming from physics, I was not aware of this keyword. I appreciate this is a difficult problem in general, but for the moment I would just be happy to understand the derivation of Przemo's result...
$endgroup$
– workandheat
Dec 3 '18 at 18:31
$begingroup$
@workandheat Now the next step is to change variables $(y,z) rightarrow (u:=(a_1,2 y+a_1,3 z)/(sqrt2a),z)$. This changes the integration region from the first quadrant to an interior of an angle $0le u < infty$ and $0le z le cdots u$. Having done that we can still do the integral over $z$ by expressing it through error functions. Finally we end up with a one dimensional integral over $u$ from one Gaussian and a product of two error functions. We do the last integral using the result from the link provided.
$endgroup$
– Przemo
Jan 23 at 11:27
add a comment |
$begingroup$
Thank you Przemo for your solution to the problem for $n=2, 3$. While I had no trouble following your derivation in 2D, I'm stuck with the derivation of your intermediate step for $n=3$. I mainly tried two approaches:
Completing the square in one variable, say $x$, leaves me with
$$int_mathbbR_+^2 mathrmdymathrmdz expleft(-frac12 fracmathrmdet,A_3mathrmdet,A_2z^2right) expleft(-frac12 fracmathrmdet, A_2a(y-m z)^2right) left[1 - mathrmerfleft(fraca_12y+a_13zsqrt2aright) right] $$
where $A_2=beginpmatrix a & a_12\ & bendpmatrix$, $A_3$ as you defined it, and $m$ is a function of the coefficients of the matrices. However, I do not know how to proceed from there: expanding the error function to do the integral in y, say, is a nightmare due to the constant term in z; I also did not find a way to do a coordinate transform à la $s=a_12y+a_13z$ or something similar.Indeed, your intermediate solution looks more like you were able to complete the square in two of the variables independently; but what happened to the cross-term? I cannot find a factorisation of the exponent that would allow me to complete two integrals over the half-line with only one variable left in the error function yielded by the integral.
Any help / hint would be greatly appreciated! Thank you in advance.
answered Dec 3 '18 at 18:09
workandheatworkandheat
587
$endgroup$
1
$begingroup$
@ workandheat, the key word here is orthant probability. If you search for this in connection with normal distribution, you get many results. But as I see even for cases with four variables people write long papers. So for more variables I doubt that there are exact results.
$endgroup$
– Karl
Dec 3 '18 at 18:23
$begingroup$
Thanks @Karl for the hint to "orthant" probability; coming from physics, I was not aware of this keyword. I appreciate this is a difficult problem in general, but for the moment I would just be happy to understand the derivation of Przemo's result...
$endgroup$
– workandheat
Dec 3 '18 at 18:31
$begingroup$
@workandheat Now the next step is to change variables $(y,z) rightarrow (u:=(a_1,2 y+a_1,3 z)/(sqrt2a),z)$. This changes the integration region from the first quadrant to an interior of an angle $0le u < infty$ and $0le z le cdots u$. Having done that we can still do the integral over $z$ by expressing it through error functions. Finally we end up with a one dimensional integral over $u$ from one Gaussian and a product of two error functions. We do the last integral using the result from the link provided.
$endgroup$
– Przemo
Jan 23 at 11:27
add a comment |
$begingroup$
Thank you Przemo for your solution to the problem for $n=2, 3$. While I had no trouble following your derivation in 2D, I'm stuck with the derivation of your intermediate step for $n=3$. I mainly tried two approaches:
Completing the square in one variable, say $x$, leaves me with
$$int_mathbbR_+^2 mathrmdymathrmdz expleft(-frac12 fracmathrmdet,A_3mathrmdet,A_2z^2right) expleft(-frac12 fracmathrmdet, A_2a(y-m z)^2right) left[1 - mathrmerfleft(fraca_12y+a_13zsqrt2aright) right] $$
where $A_2=beginpmatrix a & a_12\ & bendpmatrix$, $A_3$ as you defined it, and $m$ is a function of the coefficients of the matrices. However, I do not know how to proceed from there: expanding the error function to do the integral in y, say, is a nightmare due to the constant term in z; I also did not find a way to do a coordinate transform à la $s=a_12y+a_13z$ or something similar.Indeed, your intermediate solution looks more like you were able to complete the square in two of the variables independently; but what happened to the cross-term? I cannot find a factorisation of the exponent that would allow me to complete two integrals over the half-line with only one variable left in the error function yielded by the integral.
Any help / hint would be greatly appreciated! Thank you in advance.
answered Dec 3 '18 at 18:09
workandheatworkandheat
587
$endgroup$
Thank you Przemo for your solution to the problem for $n=2, 3$. While I had no trouble following your derivation in 2D, I'm stuck with the derivation of your intermediate step for $n=3$. I mainly tried two approaches:
Completing the square in one variable, say $x$, leaves me with
$$int_mathbbR_+^2 mathrmdymathrmdz expleft(-frac12 fracmathrmdet,A_3mathrmdet,A_2z^2right) expleft(-frac12 fracmathrmdet, A_2a(y-m z)^2right) left[1 - mathrmerfleft(fraca_12y+a_13zsqrt2aright) right] $$
where $A_2=beginpmatrix a & a_12\ & bendpmatrix$, $A_3$ as you defined it, and $m$ is a function of the coefficients of the matrices. However, I do not know how to proceed from there: expanding the error function to do the integral in y, say, is a nightmare due to the constant term in z; I also did not find a way to do a coordinate transform à la $s=a_12y+a_13z$ or something similar.Indeed, your intermediate solution looks more like you were able to complete the square in two of the variables independently; but what happened to the cross-term? I cannot find a factorisation of the exponent that would allow me to complete two integrals over the half-line with only one variable left in the error function yielded by the integral.
Any help / hint would be greatly appreciated! Thank you in advance.
answered Dec 3 '18 at 18:09
workandheatworkandheat
587
answered Dec 3 '18 at 18:09
workandheatworkandheat
587
answered Dec 3 '18 at 18:09
workandheatworkandheat
587
answered Dec 3 '18 at 18:09
workandheatworkandheat
587
587
1
$begingroup$
@ workandheat, the key word here is orthant probability. If you search for this in connection with normal distribution, you get many results. But as I see even for cases with four variables people write long papers. So for more variables I doubt that there are exact results.
$endgroup$
– Karl
Dec 3 '18 at 18:23
$begingroup$
Thanks @Karl for the hint to "orthant" probability; coming from physics, I was not aware of this keyword. I appreciate this is a difficult problem in general, but for the moment I would just be happy to understand the derivation of Przemo's result...
$endgroup$
– workandheat
Dec 3 '18 at 18:31
$begingroup$
@workandheat Now the next step is to change variables $(y,z) rightarrow (u:=(a_1,2 y+a_1,3 z)/(sqrt2a),z)$. This changes the integration region from the first quadrant to an interior of an angle $0le u < infty$ and $0le z le cdots u$. Having done that we can still do the integral over $z$ by expressing it through error functions. Finally we end up with a one dimensional integral over $u$ from one Gaussian and a product of two error functions. We do the last integral using the result from the link provided.
$endgroup$
– Przemo
Jan 23 at 11:27
add a comment |
1
$begingroup$
@ workandheat, the key word here is orthant probability. If you search for this in connection with normal distribution, you get many results. But as I see even for cases with four variables people write long papers. So for more variables I doubt that there are exact results.
$endgroup$
– Karl
Dec 3 '18 at 18:23
$begingroup$
Thanks @Karl for the hint to "orthant" probability; coming from physics, I was not aware of this keyword. I appreciate this is a difficult problem in general, but for the moment I would just be happy to understand the derivation of Przemo's result...
$endgroup$
– workandheat
Dec 3 '18 at 18:31
$begingroup$
@workandheat Now the next step is to change variables $(y,z) rightarrow (u:=(a_1,2 y+a_1,3 z)/(sqrt2a),z)$. This changes the integration region from the first quadrant to an interior of an angle $0le u < infty$ and $0le z le cdots u$. Having done that we can still do the integral over $z$ by expressing it through error functions. Finally we end up with a one dimensional integral over $u$ from one Gaussian and a product of two error functions. We do the last integral using the result from the link provided.
$endgroup$
– Przemo
Jan 23 at 11:27
1
1
$begingroup$
@ workandheat, the key word here is orthant probability. If you search for this in connection with normal distribution, you get many results. But as I see even for cases with four variables people write long papers. So for more variables I doubt that there are exact results.
$endgroup$
– Karl
Dec 3 '18 at 18:23
$begingroup$
@ workandheat, the key word here is orthant probability. If you search for this in connection with normal distribution, you get many results. But as I see even for cases with four variables people write long papers. So for more variables I doubt that there are exact results.
$endgroup$
– Karl
Dec 3 '18 at 18:23
$begingroup$
Thanks @Karl for the hint to "orthant" probability; coming from physics, I was not aware of this keyword. I appreciate this is a difficult problem in general, but for the moment I would just be happy to understand the derivation of Przemo's result...
$endgroup$
– workandheat
Dec 3 '18 at 18:31
$begingroup$
Thanks @Karl for the hint to "orthant" probability; coming from physics, I was not aware of this keyword. I appreciate this is a difficult problem in general, but for the moment I would just be happy to understand the derivation of Przemo's result...
$endgroup$
– workandheat
Dec 3 '18 at 18:31
$begingroup$
@workandheat Now the next step is to change variables $(y,z) rightarrow (u:=(a_1,2 y+a_1,3 z)/(sqrt2a),z)$. This changes the integration region from the first quadrant to an interior of an angle $0le u < infty$ and $0le z le cdots u$. Having done that we can still do the integral over $z$ by expressing it through error functions. Finally we end up with a one dimensional integral over $u$ from one Gaussian and a product of two error functions. We do the last integral using the result from the link provided.
$endgroup$
– Przemo
Jan 23 at 11:27
$begingroup$
@workandheat Now the next step is to change variables $(y,z) rightarrow (u:=(a_1,2 y+a_1,3 z)/(sqrt2a),z)$. This changes the integration region from the first quadrant to an interior of an angle $0le u < infty$ and $0le z le cdots u$. Having done that we can still do the integral over $z$ by expressing it through error functions. Finally we end up with a one dimensional integral over $u$ from one Gaussian and a product of two error functions. We do the last integral using the result from the link provided.
$endgroup$
– Przemo
Jan 23 at 11:27
add a comment |
$begingroup$
Other names for this quantity are the "multivariate Gaussian cumulative distribution", the "normalization constant of the truncated normal distribution", "non-centered orthant probabilities", ...
There appears to be a rather extensive literature on this. See for example The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting and many citations therein, like this one
Here is a paper that has closed-form expressions for the orthant probabilities for $n=4$, under different sets of assumptions for the covariance matrix.
I will update this answer as I learn more about it
answered Feb 13 at 15:18
guillefixguillefix
295111
$endgroup$
add a comment |
$begingroup$
Other names for this quantity are the "multivariate Gaussian cumulative distribution", the "normalization constant of the truncated normal distribution", "non-centered orthant probabilities", ...
There appears to be a rather extensive literature on this. See for example The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting and many citations therein, like this one
Here is a paper that has closed-form expressions for the orthant probabilities for $n=4$, under different sets of assumptions for the covariance matrix.
I will update this answer as I learn more about it
answered Feb 13 at 15:18
guillefixguillefix
295111
$endgroup$
add a comment |
$begingroup$
Other names for this quantity are the "multivariate Gaussian cumulative distribution", the "normalization constant of the truncated normal distribution", "non-centered orthant probabilities", ...
There appears to be a rather extensive literature on this. See for example The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting and many citations therein, like this one
Here is a paper that has closed-form expressions for the orthant probabilities for $n=4$, under different sets of assumptions for the covariance matrix.
I will update this answer as I learn more about it
answered Feb 13 at 15:18
guillefixguillefix
295111
$endgroup$
Other names for this quantity are the "multivariate Gaussian cumulative distribution", the "normalization constant of the truncated normal distribution", "non-centered orthant probabilities", ...
There appears to be a rather extensive literature on this. See for example The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting and many citations therein, like this one
Here is a paper that has closed-form expressions for the orthant probabilities for $n=4$, under different sets of assumptions for the covariance matrix.
I will update this answer as I learn more about it
answered Feb 13 at 15:18
guillefixguillefix
295111
edited Feb 13 at 16:53
answered Feb 13 at 15:18
guillefixguillefix
295111
answered Feb 13 at 15:18
guillefixguillefix
295111
answered Feb 13 at 15:18
guillefixguillefix
295111
295111
add a comment |
add a comment |
$begingroup$
Here we provide an answer for $n=5$ in case when the underlying matrix $bf A$ has all diagonal entries equal to unity and all cross-diagonal entries equal to $rho$ where $0le rho le 1$. We have derived the result basically in the same way as in my previous answer above, meaning by firstly bringing the quadratic form to a square in one variable and integrating over that variable and then by successively integrating over the remaining variables and reducing the dimension of the integral.
Firstly let us note that the function $mathfrak J^(1,1)$ is defined as in my previous answer above and then
let us also define the following:
beginequation
mathfrak J^(2,1)left( (a_1,a_2), b, cright):= intlimits_0^infty frace^-1/2 xi^2sqrt2 pi cdot [prodlimits_j=1^2 erf(a_j xi)] cdot T(b xi, c) dxi
endequation
This function can be always reduced to di-logarithms as shown in An integral involving a Gaussian, error functions and the Owen's T function. .
Then we have:
begineqnarray
&&P=frac2^3/2 pisqrt(1-rho)^4(1+4 rho)left(right.\
&½sqrtpi mbox$left( 2 pi arctanleft(fracrho +1sqrt1-rho ^2right)-2 arctanleft(sqrtfrac4 rho +12 rho +1right) left(3 arctanleft(fracrho +1sqrt1-rho ^2right)+pi right)+arccos(rho )
left(pi -3 arctanleft(sqrtfrac4 rho +12 rho +1right)right)right)+$\
&& 6 pi^3/2 mathfrak J^(1,1)left( fracsqrtfrac32 rho sqrt4 rho +1,fracsqrtfracrho +14 rho +1sqrt2,sqrt3right)+\
&& 6 pi^3/2 mathfrak J^(2,1)left((fracsqrt4 rho +1sqrt6 rho ,fracrho +1sqrt2 sqrt1-rho ^2),fracsqrtrho +1sqrt6 rho ,sqrt3 right)+ \
&& 6 pi^3/2 mathfrak J^(2,1)left((sqrtfrac32,sqrtfrac4 rho +1rho +1),fracsqrt6 rho sqrtrho +1,fracrho +1sqrt1-rho ^2 right)\
left.right)
endeqnarray
Below I plot the quantity in question as a function of $rho$. Note that the value $P(rho=0)= (sqrtpi/sqrt2)^5 simeq 3.09243$ as it is.
answered Mar 14 at 17:27
PrzemoPrzemo
4,40811031
$endgroup$
add a comment |
$begingroup$
Here we provide an answer for $n=5$ in case when the underlying matrix $bf A$ has all diagonal entries equal to unity and all cross-diagonal entries equal to $rho$ where $0le rho le 1$. We have derived the result basically in the same way as in my previous answer above, meaning by firstly bringing the quadratic form to a square in one variable and integrating over that variable and then by successively integrating over the remaining variables and reducing the dimension of the integral.
Firstly let us note that the function $mathfrak J^(1,1)$ is defined as in my previous answer above and then
let us also define the following:
beginequation
mathfrak J^(2,1)left( (a_1,a_2), b, cright):= intlimits_0^infty frace^-1/2 xi^2sqrt2 pi cdot [prodlimits_j=1^2 erf(a_j xi)] cdot T(b xi, c) dxi
endequation
This function can be always reduced to di-logarithms as shown in An integral involving a Gaussian, error functions and the Owen's T function. .
Then we have:
begineqnarray
&&P=frac2^3/2 pisqrt(1-rho)^4(1+4 rho)left(right.\
&½sqrtpi mbox$left( 2 pi arctanleft(fracrho +1sqrt1-rho ^2right)-2 arctanleft(sqrtfrac4 rho +12 rho +1right) left(3 arctanleft(fracrho +1sqrt1-rho ^2right)+pi right)+arccos(rho )
left(pi -3 arctanleft(sqrtfrac4 rho +12 rho +1right)right)right)+$\
&& 6 pi^3/2 mathfrak J^(1,1)left( fracsqrtfrac32 rho sqrt4 rho +1,fracsqrtfracrho +14 rho +1sqrt2,sqrt3right)+\
&& 6 pi^3/2 mathfrak J^(2,1)left((fracsqrt4 rho +1sqrt6 rho ,fracrho +1sqrt2 sqrt1-rho ^2),fracsqrtrho +1sqrt6 rho ,sqrt3 right)+ \
&& 6 pi^3/2 mathfrak J^(2,1)left((sqrtfrac32,sqrtfrac4 rho +1rho +1),fracsqrt6 rho sqrtrho +1,fracrho +1sqrt1-rho ^2 right)\
left.right)
endeqnarray
Below I plot the quantity in question as a function of $rho$. Note that the value $P(rho=0)= (sqrtpi/sqrt2)^5 simeq 3.09243$ as it is.
answered Mar 14 at 17:27
PrzemoPrzemo
4,40811031
$endgroup$
add a comment |
$begingroup$
Here we provide an answer for $n=5$ in case when the underlying matrix $bf A$ has all diagonal entries equal to unity and all cross-diagonal entries equal to $rho$ where $0le rho le 1$. We have derived the result basically in the same way as in my previous answer above, meaning by firstly bringing the quadratic form to a square in one variable and integrating over that variable and then by successively integrating over the remaining variables and reducing the dimension of the integral.
Firstly let us note that the function $mathfrak J^(1,1)$ is defined as in my previous answer above and then
let us also define the following:
beginequation
mathfrak J^(2,1)left( (a_1,a_2), b, cright):= intlimits_0^infty frace^-1/2 xi^2sqrt2 pi cdot [prodlimits_j=1^2 erf(a_j xi)] cdot T(b xi, c) dxi
endequation
This function can be always reduced to di-logarithms as shown in An integral involving a Gaussian, error functions and the Owen's T function. .
Then we have:
begineqnarray
&&P=frac2^3/2 pisqrt(1-rho)^4(1+4 rho)left(right.\
&½sqrtpi mbox$left( 2 pi arctanleft(fracrho +1sqrt1-rho ^2right)-2 arctanleft(sqrtfrac4 rho +12 rho +1right) left(3 arctanleft(fracrho +1sqrt1-rho ^2right)+pi right)+arccos(rho )
left(pi -3 arctanleft(sqrtfrac4 rho +12 rho +1right)right)right)+$\
&& 6 pi^3/2 mathfrak J^(1,1)left( fracsqrtfrac32 rho sqrt4 rho +1,fracsqrtfracrho +14 rho +1sqrt2,sqrt3right)+\
&& 6 pi^3/2 mathfrak J^(2,1)left((fracsqrt4 rho +1sqrt6 rho ,fracrho +1sqrt2 sqrt1-rho ^2),fracsqrtrho +1sqrt6 rho ,sqrt3 right)+ \
&& 6 pi^3/2 mathfrak J^(2,1)left((sqrtfrac32,sqrtfrac4 rho +1rho +1),fracsqrt6 rho sqrtrho +1,fracrho +1sqrt1-rho ^2 right)\
left.right)
endeqnarray
Below I plot the quantity in question as a function of $rho$. Note that the value $P(rho=0)= (sqrtpi/sqrt2)^5 simeq 3.09243$ as it is.
answered Mar 14 at 17:27
PrzemoPrzemo
4,40811031
$endgroup$
Here we provide an answer for $n=5$ in case when the underlying matrix $bf A$ has all diagonal entries equal to unity and all cross-diagonal entries equal to $rho$ where $0le rho le 1$. We have derived the result basically in the same way as in my previous answer above, meaning by firstly bringing the quadratic form to a square in one variable and integrating over that variable and then by successively integrating over the remaining variables and reducing the dimension of the integral.
Firstly let us note that the function $mathfrak J^(1,1)$ is defined as in my previous answer above and then
let us also define the following:
beginequation
mathfrak J^(2,1)left( (a_1,a_2), b, cright):= intlimits_0^infty frace^-1/2 xi^2sqrt2 pi cdot [prodlimits_j=1^2 erf(a_j xi)] cdot T(b xi, c) dxi
endequation
This function can be always reduced to di-logarithms as shown in An integral involving a Gaussian, error functions and the Owen's T function. .
Then we have:
begineqnarray
&&P=frac2^3/2 pisqrt(1-rho)^4(1+4 rho)left(right.\
&½sqrtpi mbox$left( 2 pi arctanleft(fracrho +1sqrt1-rho ^2right)-2 arctanleft(sqrtfrac4 rho +12 rho +1right) left(3 arctanleft(fracrho +1sqrt1-rho ^2right)+pi right)+arccos(rho )
left(pi -3 arctanleft(sqrtfrac4 rho +12 rho +1right)right)right)+$\
&& 6 pi^3/2 mathfrak J^(1,1)left( fracsqrtfrac32 rho sqrt4 rho +1,fracsqrtfracrho +14 rho +1sqrt2,sqrt3right)+\
&& 6 pi^3/2 mathfrak J^(2,1)left((fracsqrt4 rho +1sqrt6 rho ,fracrho +1sqrt2 sqrt1-rho ^2),fracsqrtrho +1sqrt6 rho ,sqrt3 right)+ \
&& 6 pi^3/2 mathfrak J^(2,1)left((sqrtfrac32,sqrtfrac4 rho +1rho +1),fracsqrt6 rho sqrtrho +1,fracrho +1sqrt1-rho ^2 right)\
left.right)
endeqnarray
Below I plot the quantity in question as a function of $rho$. Note that the value $P(rho=0)= (sqrtpi/sqrt2)^5 simeq 3.09243$ as it is.
answered Mar 14 at 17:27
PrzemoPrzemo
4,40811031
answered Mar 14 at 17:27
PrzemoPrzemo
4,40811031
answered Mar 14 at 17:27
PrzemoPrzemo
4,40811031
answered Mar 14 at 17:27
PrzemoPrzemo
4,40811031
4,40811031
add a comment |
add a comment |
Thanks for contributing an answer to Mathematics Stack Exchange!
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
Use MathJax to format equations. MathJax reference.
To learn more, see our tips on writing great answers.
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f869502%2fmultivariate-gaussian-integral-over-positive-reals%23new-answer', 'question_page');
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
$begingroup$
Your notation, $mathbf xgeq0$ doesn't make sense when $mathbfxinmathbbR^n$, $ngeq 2$.
$endgroup$
– nullgeppetto
Jul 17 '14 at 14:21
$begingroup$
I think that you should use the solution provided in the link you give about my older question.
$endgroup$
– nullgeppetto
Jul 17 '14 at 14:23
$begingroup$
With $x geq 0$ I meant $x_i geq 0$ $forall i$. Is this bad notation? Your question is different from mine in that your integral can be reduced to n-1 gaussian integrals and one "half gaussian" one, which can be solved by the error function. I am not aware of an n-dimensional error function though...
$endgroup$
– le_m
Jul 17 '14 at 14:32
$begingroup$
I encountered similar problem recently. For the case of all covariances being the same, the problem is exactly solvable. If this is the one you are interested in, let me know.
$endgroup$
– Sungmin
Nov 6 '15 at 14:55
$begingroup$
@Sungmin: Would you mind showing me the solution in this case ? Honestly I truly doubt that any results exist when $n > 4$. For $n=4$ di-logarithms enter the result and in case of higher dimensions higher order poly-loarithms will be needed.
$endgroup$
– Przemo
Mar 12 at 18:35