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]?
 A: 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.
A: 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()

