# Calculating a conditional expectation

Suppose that $$X$$ follows a beta distribution with parameters $$\alpha > 0$$, $$\beta > 0$$. I want to calculate $$\frac{\int_c^1 x f(x) dx}{P(X \geq c)} - \frac{\int_0^c x f(x) dx}{P(X \leq c)}$$

for some constant $$c \in (0, 1)$$.

I have tried the following code. While it works in the $$\alpha = \beta = 1$$ case, it seems won't run (at least in a reasonable amount of time) in other cases. What am I doing wrong? The code is below:

less[c_, \[Alpha]_, \[Beta]_] :=
Integrate[x PDF[BetaDistribution[\[Alpha], \[Beta]], x], {x, 0, c}]/
CDF[BetaDistribution[\[Alpha], \[Beta]], c]

more[c_, \[Alpha]_, \[Beta]_] :=
Integrate[
x PDF[BetaDistribution[\[Alpha], \[Beta]], x], {x, c, 1}]/(1 -
CDF[BetaDistribution[\[Alpha], \[Beta]], c])

diff[c_, \[Alpha]_, \[Beta]_] :=
more[c, \[Alpha], \[Beta]] - less[c, \[Alpha], \[Beta]]

Plot[diff[c, 0.5, 0.5], {c, 0, 1}]

• I think that you can calculate all the integrals in closed form; that might help with performance. Jan 30 at 15:23

## 1 Answer

\$Version

(* "13.2.0 for Mac OS X x86 (64-bit) (November 18, 2022)" *)

Clear["Global*"]

dist = BetaDistribution[α, β];

dpa = DistributionParameterAssumptions[dist]

(* α > 0 && β > 0 *)

DistributionDomain[dist]

(* Interval[{0, 1}] *)


Integrals are generally faster and simpler when given the appropriate assumptions.

less[c_, α_, β_] = Assuming[dpa && 0 < c < 1,
Integrate[x PDF[dist, x], {x, 0, c}]/CDF[dist, c] // Simplify]

(* Beta[c, 1 + α, β]/(
Beta[α, β] BetaRegularized[c, α, β]) *)

more[c_, α_, β_] = Assuming[dpa && 0 < c < 1,
Integrate[x PDF[dist, x], {x, c, 1}]/(1 - CDF[dist, c]) // Simplify]

(* (-Gamma[1 + α] Gamma[β] +
Beta[c, 1 + α, β] Gamma[
1 + α + β])/(Beta[α, β] (-1 +
BetaRegularized[c, α, β]) Gamma[
1 + α + β]) *)


EDIT: Use FunctionExpand

diff[c_, α_, β_] = Assuming[dpa && 0 < c < 1,
more[c, α, β] - less[c, α, β] //
FunctionExpand // FullSimplify]

(* ((1 - c)^β c^α Gamma[α] \
Gamma[β])/((α + β) Beta[
c, α, β] Gamma[α] Gamma[β] -
Beta[c, α, β]^2 Gamma[1 + α + β]) *)

Plot[diff[c, 0.5, 0.5], {c, 0, 1}]


EDIT: Or more simply,

p[c_, α_, β_] = Assuming[0 < c < 1,
Expectation[x \[Conditioned] x >= c, x \[Distributed] dist] -
Expectation[x \[Conditioned] x < c, x \[Distributed] dist] //
FullSimplify]

(* ((1 - c)^β c^α Gamma[α] Gamma[β])/((α + \
β) Beta[c, α, β] Gamma[α] Gamma[β] -
Beta[c, α, β]^2 Gamma[1 + α + β]) *)

diff[c, α, β] === p[c, α, β]]

(* True *)
`