# How to help MMA to simplify integrands?

MMA sometimes has to grind for LONG time on easy integrals:

gau[x_, m_, var_] := (E^(-(x - m)^2/(2 var)));
f[x_] := gau[x, x0, a];
p[x_] := gau[x, x1, b];
q[x_] := gau[x, x2, c];
integrand = Simplify[f[x]*q[x]/p[x]];
Integrate[integrand, {x, -Infinity, Infinity}]


With a little algebra, this just comes the definite integral of a gaussian. I know I could to all the exponent algebra before the Integrate, but I can't always anticipate the situation. Any tips?

Just to be clear, my code above WORKS, but it takes so long for such a simple substitution and look-up that I suspect something could be done to narrow MMA's search.

• The integral took 14 seconds on my laptop (V10). Is that really a problem? – Sjoerd C. de Vries Nov 9 '14 at 21:27
• @SjoerdC.deVries Yes, because I have to do hundreds of them. – Jerry Guern Nov 9 '14 at 21:29
• Integration is, as you will be quite aware of, not a straightforward process. There is not a single path that is guaranteed to lead to the solution and many different tricks may have to be used before a solution is reached. I note that the indefinite integral in this case is found almost instantaneously and I assume that Mathematica spends most time in determining appropriate conditions for the constants. And then, a few hundred integrals taking the timing I got should not be really a problem. Looks like they would have been calculated in the time you spent here on this question ;-) – Sjoerd C. de Vries Nov 9 '14 at 22:09
• "I can't always anticipate the situation" -- I'm not sure what you mean. Do you mean any random function is possible? Are they always gaussian, perhaps times a constant? Somewhere between these two extremes? One certainly can do Gaussians in less than 0.001 sec.... – Michael E2 Jan 9 '15 at 2:07
• @MichaelE2 f,p,&q will always be normalized pdfs, and fq/p will always be integrable. f will usually be Gaussian, not always, and q&p will sometimes be Gaussian mixtures sometimes with one term each. So I could write code to simplify fq/p into a single gaussian and get the integral analytically, but I won't know in advance when f,p&q are all gaussians. – Jerry Guern Jan 9 '15 at 7:08

No Problem in Version 8.

\$Version

(* Out[114]= "8.0 for Microsoft Windows (64-bit) (October 7, 2011)" *)

gau[x_, m_, var_] := (E^(-(x - m)^2/(2 var)))

f[x_] := gau[x, x0, a]

p[x_] := gau[x, x1, b]

q[x_] := gau[x, x2, c]

integrand = Simplify[f[x]*q[x]/p[x]]

(* Out[112]= E^(1/2 (-((x - x0)^2/a) + (x - x1)^2/b - (x - x2)^2/c)) *)

Integrate[integrand, {x, -Infinity, Infinity}]

(* Out[108]= ConditionalExpression[(
E^((c (x0 - x1)^2 - b (x0 - x2)^2 + a (x1 - x2)^2)/(2 (a (b - c) + b c))) Sqrt[2 \[Pi]])/Sqrt[1/a - 1/b + 1/c], Re[1/a - 1/b + 1/c] > 0] *)


Regards,
Wolfgang

It's hard to know what exactly will help, given the possible variations in the functions. One suggestion is to use GenerateConditions -> False, so that time is not spent on figuring out the complicated conditions for which the value of the integral is valid, as Sjoerd surmises; time will still be spent trying to figure out the limits at infinity (at least, I'm pretty sure that is a complicated calculation if you're going to worry about convergence). Another idea is to use Expectation, which is fast on certain expressions when x is normally distributed. The function int below should work and be fast when all three functions are Gaussian, as in the example. With some more work, one could split a Gaussian mixture pdf into its Gaussian components, integrate, and add. One might have to test what happens when f q / p is not reducible to Gaussians. Expectation might default to Integration such cases. Also, for some inputs the NormalDistribution constructed by int might be invalid.

Timings:

Integrate[integrand, {x, -Infinity, Infinity}]; // AbsoluteTiming
Integrate[integrand, {x, -Infinity, Infinity},
GenerateConditions -> False]; // AbsoluteTiming
int[f[x], p[x], q[x], x]; // AbsoluteTiming
(*
{20.370266, Null}
{5.157397, Null}
{0.001574, Null}
*)


Function:

int::nnormal = "1 resulted in an invalid normal distribution. Using Integrate.";
int[f_, p_, q_, x_] := Module[{factor, exponent, a, b, c, integrand},
integrand = Simplify[f*q/p];
If[MatchQ[integrand, Power[E, _]],
factor = 1,
If[MatchQ[integrand, _Times],
factor = Times @@ Cases[integrand, Except[Power[E, _]]],
factor = integrand (* fails *)
]];
(*If[InternalDependsOnQ[factor, x],
Return[Integrate[integrand, {x, -Infinity, Infinity},
GenerateConditions -> False, Assumptions -> IntegrategetAllVariables[integrand, x] ∈ Reals]]
];*)
exponent = Simplify[integrand/factor] /. Power[E, z_] :> z;
{c, b, a} = PadRight[CoefficientList[exponent, x], 3];
If[a =!= 0,
Sqrt[π]/Sqrt[-a] Exp[a (c/a - (b/(2 a))^2)] *
Expectation[factor, x \[Distributed] NormalDistribution[-b/(2 a), Sqrt[-1/(2 a)]]],
Message[int::nnormal, integrand];
Integrate[integrand, {x, -Infinity, Infinity},
GenerateConditions -> False, Assumptions -> IntegrategetAllVariables[integrand, x] ∈ Reals]
]
];


Uncomment the test InternalDependsOnQ[factor, x] if you want the computation to default to Integrate when the factor is not a constant. In that case,

Expectation[factor, x \[Distributed] NormalDistribution[-b/(2 a), Sqrt[-1/(2 a)]]]


evaluates to factor, and the code could be shortened.

For alternatives to IntegrategetAllVariables[integrand, x], see Extracting variables from an expression and Get all variables in an expression with Variables[]. Sometimes it helps to include x ∈ Reals in assumptions, which is not included in the following:

IntegrategetAllVariables[integrand, x]
(*  {a, b, c, x0, x1, x2}  *)