0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ The integral took 14 seconds on my laptop (V10). Is that really a problem? $\endgroup$ – Sjoerd C. de Vries Nov 9 '14 at 21:27
  • $\begingroup$ @SjoerdC.deVries Yes, because I have to do hundreds of them. $\endgroup$ – Jerry Guern Nov 9 '14 at 21:29
  • 3
    $\begingroup$ 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 ;-) $\endgroup$ – Sjoerd C. de Vries Nov 9 '14 at 22:09
  • $\begingroup$ "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.... $\endgroup$ – Michael E2 Jan 9 '15 at 2:07
  • $\begingroup$ @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. $\endgroup$ – Jerry Guern Jan 9 '15 at 7:08
1
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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}  *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.