4
$\begingroup$

I'm using Expectation to calculate the Gaussian integral of a user-supplied function. The following works well and fast (< 1 second):

a[xi_, xj_] := E^(-1/2*(xi - xj)^2/σa^2);

Expectation[a[x[i], x[j]], x[i] \[Distributed] NormalDistribution[xav[i], σ[i]]]

enter image description here

If I try it on an indefinite Sum, it takes 165 seconds to ultimately fail:

Expectation[
  Sum[a[x[i], x[j]], {j, n}, Method -> "Procedural"],
  x[i] \[Distributed] NormalDistribution[xav[i], σ[i]]
]

enter image description here

Manually switching the order to a Sum of the Expectation works great (again, < 1 sec):

Sum[
  Expectation[a[x[i], x[j]], x[i] \[Distributed] NormalDistribution[xav[i], σ[i]]],
  {j, n}, Method -> "Procedural"
]

enter image description here

However, all this is happening inside another function that takes arbitrary input (including the Sum), so I want to switch the Sum and Expectation automatically. The following replacement rule does the trick, but is very slow (165 seconds):

Expectation[
  Sum[a[x[i], x[j]], {j, n}, Method -> "Procedural"],
  x[i] \[Distributed] NormalDistribution[xav[i], σ[i]]
] /. Expectation[Sum[func_, range_, opts___], dist_] :> Sum[Expectation[func, dist], range, opts]

enter image description here

Is there a better (faster, more elegant) way to achieve the same output?

Improve this question
$\endgroup$
3
  • 4
    $\begingroup$ I'm a bit confused by the notation. If x[i] is a random variable, then isn't x[j] also a random variable? If so, the distribution of x[j] seems to be ignored. Also, because the sum over j includes i, the case ofa[x[i],x[i]] occurs. The value of Expectation[a[x[i], x[i]], x[i] \[Distributed] NormalDistribution[xav[i], \[Sigma][i]]] becomes 1 which is not the general form that you show. Am I missing something? $\endgroup$
    – JimB
    Mar 26 at 4:25
  • $\begingroup$ Thanks for the comments, which helped me realize that I set up my toy example wrong (should be the x[j] that are distributed not the x[i]). The problem still exists, but I'll need to put the Expectations inside the Sum earlier, so maybe not relevant to me anymore. $\endgroup$
    – Chris K
    Mar 26 at 17:40
  • $\begingroup$ The problem comes from ecological models with intraspecific trait variation. The x[i] and x[j] are the traits of interacting individuals, who compete with competition kernel a[x[i], x[j]]. Assuming each species' traits are normally distributed, we have to integrate over their distributions to compute competition. See this paper for more info. $\endgroup$
    – Chris K
    Mar 26 at 17:40

3 Answers 3

Reset to default
4
$\begingroup$

This is an extended comment.

If all of the x[i] are random variables with NormalDistribution[xav[i], σ[i]], I wonder if what you really want is the following:

a[xi_, xj_] := E^(-1/2*(xi - xj)^2/\[Sigma]a^2);

If i != j, then we have

FullSimplify[Expectation[a[x[i], x[j]], {x[i] \[Distributed] NormalDistribution[xav[i], σ[i]],
  x[j] \[Distributed] NormalDistribution[xav[j], σ[j]]}], 
  Assumptions -> {σa > 0, σ[i] > 0, σ[j] > 0}]

Expectation for i and j

But if i=j, then one has

FullSimplify[Expectation[a[x[i], x[i]], x[i] \[Distributed] NormalDistribution[xav[i], σ[i]]], 
  Assumptions -> {σa > 0, σ[i] > 0}]

(* 1 *)

So you could just define the expectation of a[xi_, xj_] as

ea[i_, j_] := Piecewise[{{1,i==j}},(E^(-((xav[i] - xav[j])^2/(2 (σa^2 + σ[i]^2 + σ[j]^2)))) σa)/Sqrt[σa^2 + σ[i]^2 + σ[j]^2]]

and

Sum[ea[i, j], {j, n}]

Sum of ea values

Improve this answer
$\endgroup$
3
  • $\begingroup$ Very insightful, thanks! I shouldn't post questions at 5pm on Friday afternoon ;) $\endgroup$
    – Chris K
    Mar 26 at 17:42
  • $\begingroup$ I know the feeling well (and it can be even worse at 10pm especially for answers). I wonder if you should be considering correlated normals as opposed to independent normals. $\endgroup$
    – JimB
    Mar 26 at 18:40
  • $\begingroup$ Good thought, but I think independent normals is appropriate, since the distributions represent separate species' phenotypic distributions. $\endgroup$
    – Chris K
    Mar 27 at 16:30
6
$\begingroup$

Include assumptions

Clear["Global`*"]

a[xi_, xj_] := E^(-1/2*(xi - xj)^2/σa^2);

Assuming[
  {σa > 0, Element[n, PositiveIntegers]},
  Expectation[Sum[a[x[i], x[j]], {j, n}], 
   x[i] \[Distributed] 
    NormalDistribution[xav[i], σ[i]]]] // AbsoluteTiming

(* {10.6525, (E^(-((x[j] - xav[i])^2/(2 (σa^2 + σ[i]^2))))
   n σa)/Sqrt[σa^2 + σ[i]^2]} *)
Improve this answer
$\endgroup$
4
  • $\begingroup$ Interesting! Although 10 seconds is still slower than my manual switch time of 1, and Method -> "Procedural" really hurts (500 s!) $\endgroup$
    – Chris K
    Mar 26 at 3:20
  • $\begingroup$ Oops, also not the right answer, since it seems to assume that all x[j] are the same! $\endgroup$
    – Chris K
    Mar 26 at 3:37
  • 1
    $\begingroup$ You only showed one of the x[j] as random (specifically x[i]). You need to reformulate the problem if all of the x[j] are random. $\endgroup$
    – Bob Hanlon
    Mar 26 at 3:53
  • $\begingroup$ Sorry, I didn't notice your last comment before I wrote my comment on the question. $\endgroup$
    – JimB
    Mar 26 at 4:28
4
$\begingroup$ 
sum /: Expectation[sum[a_, b__], c___] := sum[Expectation[a, c], b]

ReplaceAll[sum -> Sum]@
  Block[{Sum = sum}, Expectation[Sum[a[x[i], x[j]], {j, n}, Method -> "Procedural"], 
    x[i] \[Distributed] NormalDistribution[xav[i], σ[i]]]] // AbsoluteTiming

enter image description here

$Version

"11.3.0 for Microsoft Windows (64-bit) (March 7, 2018)"
Improve this answer
$\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.