# Characteristic Function of 3D Gaussian with Boole

Context

Given the Probability Distribution

pdf=ProbabilityDistribution[(1/(8 π))675 Sqrt[5] *
(x1-x2) (x1-x3) (x2-x3)
Exp[1/2 (-x1 (6 x1-(3 x2)/2-(3 x3)/2)-x2 (-((3 x1)/2)+6 x2-(3 x3)/2)-
x3 (-((3 x1)/2)-(3 x2)/2+6 x3))],
{x1,-∞,∞},{x2,-∞,x1},{x3,-∞,x2}]


Note the Boundary on x2 and x3.

Question

I would like to compute the Characteristic Function of this PDF.

Attempt

 CF=CharacteristicFunction[pdf, {y1,y2,y3}]


takes forever.

I have tried rectifying the boundary as follows:

rs = Solve[{(x3 - x2)/2 == a, (x2 - x1)/2 == b, x3 == c}, {x1, x2, x3}][[1]];
pdf2 = Exp[I x1 y1 + I x2 y2 + I x3 y3 ] PDF[pdf,{x1,x2,x3}] /. rs // Simplify


(* (-675*Sqrt[5]ab*(a + b)* E^(-18*a^2 - 12*b^2 - (9*c^2)/2 + b*(6*c - (2*I)*y1) - 2*a*(9*b - 6*c + I*(y1 + y2)) + I * c*(y1 + y2 + y3)))/π *)

But this does not seem to converge either:

Integrate[pdf2 , {c, -∞, ∞}, {b, -∞,0}, {a, -∞, 0}]


runs forever.

I can integrate over 2 out of the 3 variables.

Integrate[pdf2 , {b, -∞, 0}, {c, -∞, ∞}]


but last integral also runs forever.

Any guidance would be welcome! :-)

Note: The purpose of knowing this Characteristic Function is to compute symbolically the matrix of scalar products involved in this question, but for 3D curvature.

• You could ask the genius integrators over at math.stackexchange.com Jan 21 '19 at 9:03
• Don't forget the factor of 4 coming from the Jacobian of your coordinate transformation. As pdf2 stands now, its integral is equal to $1/4$ instead of 1. You can find the Jacobian with D[{x1, x2, x3} /. rs, {{a, b, c}}] // Det // Abs. Jan 21 '19 at 9:08
• Thanks for your comments Jan 21 '19 at 17:55

## 1 Answer

Following the suggestion of this answer, the Characteristic function can be partially computed after the following change of variable on the PDF:

rs = Solve[{(x3 - x2)/2 == a, (x2 - x1)/2 == b, x3 == c}, {x1, x2, x3}][[1]]


as

Pdf2 = Exp[I x1 y1 + I x2 y2 + I x3 y3 ] Pdf1 /. rs // Simplify


(* -((675 Sqrt[5] a b (a+b) exp(-18 a^2-2 a (9 b-6 c+I (y1+y2))-12 b^2+b (6 c-2 I y1)-(9 c^2)/2+I c (y1+y2+y3)))/π) *)

Then the first two integrals can be done.

CF2 = 4 Integrate[Pdf2 , {b, -Infinity, 0}, {c, -Infinity, Infinity}]


The good news is that if you only want the moments (as it happens) they can still be computed exactly as follows

Clear[mom];
mom[p1_, p2_, p3_] :=
mom[p1, p2, p3] =
(1/I)^(p1 + p2 + p3) D[CF2, {y1, p1}, {y2, p2}, {y3, p3}] /.
Thread[{y1, y2, y3} -> 0] //
FullSimplify[#, Assumptions -> a < 0] & //
Integrate[#, {a, -Infinity, 0}] &


FYI they involve quantities such as

4/375 (283 Sqrt[10/π]-32 Sqrt[30 π])