How to help MMA to simplify integrands?

Question

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.

Answer 1

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

Answer 2

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[Internal`DependsOnQ[factor, x],
    Return[Integrate[integrand, {x, -Infinity, Infinity},
     GenerateConditions -> False, Assumptions -> Integrate`getAllVariables[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 -> Integrate`getAllVariables[integrand, x] ∈ Reals]
    ]
   ];

Uncomment the test Internal`DependsOnQ[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 Integrate`getAllVariables[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:

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