# How to calculate a partial expected value of beta distribution (mean of a truncated beta)?

Given a Beta Distribution with a=2, b=3, we can find an expected value (mean) for the interval [0, 1] = a/(a+b) = 2/5 = 0.4 and median = (a - 1/3)/(a+b-2/3) = 0.39, which are close.

I am looking for a solution in python. I can use scipy.stats.beta to calculate median for the interval [0, 0.4] with percent point function (inverse of cdf — percentiles):

beta.ppf(0.4/2,a,b) = 0.2504


Since for this beta distribution, the overall mean and median are close (0.4 and 0.39 respectively), I use the median for the interval [0, 0.4] to estimate the expected values (mean) for the interval [0, 0.4].

Is there any way to calculate expected values (mean) for the interval [0, 0.4]?

• As a general rule, you should not correct a calculation error in your question that is explained in an answer, because you "break" the answer - the correction in the answer no longer makes sense. (On the other hand making a correction offered in comments is a different matter.) --- I try to will adjust my answer to compensate but such errors are usually best left as is. Commented Sep 7, 2017 at 2:06
• Sorry about that. I didn't know such rule. Commented Sep 7, 2017 at 2:24

Note that the formula you have near the top there for the beta median ($\frac{\alpha-\frac13}{\alpha+\beta-\frac23}$) is approximate. You should be able to compute an effectively "exact" numerical median with the inverse cdf (quantile function) of the beta distribution in Python (for a $\text{beta}(2,3)$ I get a median of around $0.3857$ while that approximate formula gives $0.3846$).

This mean of a truncated distribution is pretty straightforward with a beta. For a positive random variable we have

$E(X|X<k) = \int_0^k x\,f(x)\, dx / \int_0^k f(x)\, dx$

where in this case $f$ is the density of a beta with parameters $\alpha$ and $\beta$ (which I'll now write as $f(x;\alpha,\beta)$):

$f(x;\alpha,\beta)=\frac{1}{B(\alpha,\beta)} x^{\alpha-1} (1-x)^{\beta-1}\,,\:0<x<1,\alpha,\beta>0$

Hence $x\,f(x) = \frac{B(\alpha+1,\beta)}{B(\alpha,\beta)} f(x;\alpha+1,\beta)=\frac{\alpha}{\alpha+\beta} f(x;\alpha+1,\beta)$

So $E(X|X<k) = \frac{\alpha}{\alpha+\beta} \int_0^k f(x;\alpha+1,\beta)\, dx / \int_0^k f(x;\alpha,\beta)\,dx$

Now the two integrals are just beta CDFs which you have available in Python already.

With $\alpha=2,\beta=3,k=0.4$ we get $E(X|X<0.4)\approx 0.24195$. This is consistent with simulation ($10^6$ simulations gave $\approx 0.24194$).

For the median, I get $F^{-1}(\frac12 F(0.4;2,3);2,3)\approx 0.25040$, which is again consistent with simulation ($10^6$ simulations gave $\approx 0.25038$).

The two are pretty close in this case but that's not a general result; they may sometimes be more substantially different.

• Thanks for your detailed explanation. I should have asked here weeks ago! Commented Sep 7, 2017 at 1:40

For people who just want a scipy one-liner, it is also possible to use the expect function:

from scipy.stats import beta as B

start = 0
end = 0.4
a=2
b=3

print(B(a, b, loc=0, scale=1).expect(lb=start,ub=end,conditional=True))


0.24195121951219523

This has the advantage of being (slightly) more readable, and more general

However, it seems significantly slower (~4X slower) than the answer of @Glen_b

Code used to generate the figure :

from scipy.stats import beta as B
from datetime import datetime
import time
import matplotlib.pyplot as plt

N_tests = 10000
start = 0
end = 0.4
a = 2
b = 3

for i in range(N_tests):
t1 = datetime.now()
value = B(a=a, b=b, loc=0, scale=1).expect(lb=start,ub=end,conditional=True)
t2 = datetime.now()
list_times1.append((t2 - t1).total_seconds())

for i in range(N_tests):
t1 = datetime.now()
value = (a/(a+b)) * B(a+1,b).cdf(end) / B(a,b).cdf(end)
t2 = datetime.now()
list_times2.append((t2 - t1).total_seconds())

_ = plt.hist(list_times2,bins=100,label="scipy cdf",density=True)
_ = plt.hist(list_times1,bins=100,label="scipy expect",density=True)
plt.xlabel("Time taken")
plt.ylabel("Number of occurences")
plt.title("Distribution of the execution times for 10000 computations")
plt.legend()

plt.savefig("expected_partial_beta.png")

plt.show()