Suppose $\theta \sim \beta(\alpha_1,\alpha_2)$ and we believe that $E[\theta] = m$ and $P(l < \theta < u) = 0.95$. Write a program that can solve for $\alpha_1$ and $\alpha_2$ in terms of $m$, $l$, and $u$.
The mean for the beta distribution is well-known, so we have that
\begin{equation} E[\theta] = m = \frac{\alpha_1}{\alpha_1 + \alpha_2}. \end{equation}Thus, this implies that
\begin{equation} \alpha_2 = \alpha_1\frac{1-m}{m}. \end{equation}Now, as $\alpha_1$ increases the strength of our prior increases, so $P(l < \theta < u)$ is a monotonic function of $\alpha_1$. And so, we can solve for $\alpha_1$, and thereby, $\alpha_2$ with a binary search.
from scipy import stats
m = 0.15
l = 0.05
u = 0.3
alpha_1_lower = 0.01
alpha_1_upper = 1000
alpha_1 = (alpha_1_lower + alpha_1_upper)/2
tol = 1e-15
while alpha_1_upper - alpha_1_lower >= 1e-9:
alpha_2 = alpha_1*(1-m)/m
density = stats.beta.cdf(u, a = alpha_1, b = alpha_2) - stats.beta.cdf(l, a = alpha_1, b = alpha_2)
if density > 0.95:
alpha_1_upper = alpha_1
else:
alpha_1_lower = alpha_1
alpha_1 = (alpha_1_lower + alpha_1_upper)/2
print(alpha_1)
print(alpha_2)
4.506062413888118 25.534353681276208
The problem asks for $\alpha_1$ and $\alpha_2$ if $m = 0.15$, $l = 0.05$, and $u = 0.3$. We find that
\begin{align} \alpha_1 &\approx 4.506062413888118 \\ \alpha_2 &\approx 25.534353681276208. \end{align}