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.

The integral took 14 seconds on my laptop (V10). Is that really a problem? — Sjoerd C. de Vries, Nov 09 '14 at 21:27
@SjoerdC.deVries Yes, because I have to do hundreds of them. — Jerry Guern, Nov 09 '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 09 '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 09 '15 at 02: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 09 '15 at 07:08
It seems that Mma makes some algebra in the initial stage of the process of calculating integrals. Am I not right? — Alexei Boulbitch, Jan 09 '15 at 09:50
@AlexeiBoulbitch It must, but if it was doing the algebra in this problem to combine the gaussians into one gaussians, this prob would get done in a fraction of a second instead of several seconds. — Jerry Guern, Jan 09 '15 at 11:58
It seems that the time taken by the operation is related first of all to the number of the parameters like a, b and c and is spent for checking relations between them. This Timing[Integrate[integrand, {x, -Infinity, Infinity}]]returns 13.5 s, the same with the assumption: Assumptions -> {a > 0, b > 0, c > 0, b > a, b > c}gives 7.2 s, while if instead I replace the variables within the integrand: Timing[Integrate[ integrand /. {a -> 1, b -> 2, c -> 1}, {x, -Infinity, Infinity}]]it yields 1.3 s. — Alexei Boulbitch, Jan 09 '15 at 13:54

score 1 · Answer 1 · answered Nov 09 '14 at 21:18

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

score 1 · Answer 2 · edited Apr 13 '17 at 12:55

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}  *)

How to help MMA to simplify integrands?

2 Answers2

Linked