# Expectation of a Sum

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]]]


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]]
]


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"
]


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]


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

• 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?
– JimB
Mar 26 at 4:25
• 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. Mar 26 at 17:40
• 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. Mar 26 at 17:40

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}]


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}]


• Very insightful, thanks! I shouldn't post questions at 5pm on Friday afternoon ;) Mar 26 at 17:42
• 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.
– JimB
Mar 26 at 18:40
• Good thought, but I think independent normals is appropriate, since the distributions represent separate species' phenotypic distributions. Mar 27 at 16:30

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

• Interesting! Although 10 seconds is still slower than my manual switch time of 1, and Method -> "Procedural" really hurts (500 s!) Mar 26 at 3:20
• Oops, also not the right answer, since it seems to assume that all x[j] are the same! Mar 26 at 3:37
• 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. Mar 26 at 3:53
• Sorry, I didn't notice your last comment before I wrote my comment on the question.
– JimB
Mar 26 at 4:28
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


\$Version

"11.3.0 for Microsoft Windows (64-bit) (March 7, 2018)"
`