alpha and beta using mean and variance of beta distribution

It is well known that for beta(alpha, beta), the mean is alpha/(alpha+beta) and the variance is (alpha*beta)/((alpha+beta)^2*(alpha+beta+1)).

https://en.wikipedia.org/wiki/Beta_distribution

Then it is easy to reverse the model to estimate alpha and beta when mean and variance are known:

alpha/p = alpha+beta
beta = alpha/p - alpha = alpha * (1/p -1)

v = (alpha*(alpha/p - alpha)) / ((alpha + alpha/p - alpha)^2 * (alpha + alpha/p - alpha +1))
v = (alpha*(alpha * (1/p -1))) / ((alpha + alpha/p - alpha)^2 * (alpha + alpha/p - alpha +1))
v = (alpha^2 * (1/p -1)) / ((alpha^2 (1/p^2)) * (alpha/p + 1))
v = (alpha * (1/p -1)) / ((alpha (1/p^2)) * (alpha/p + 1))
v = (1/p -1) / (((1/p^2)) * (alpha/p + 1))
1/v = ((1/p^2) * (alpha/p + 1)) / (1/p -1)
(1/p -1) / v = (1/p^2) * (alpha/p + 1)
(1/p -1) / ( v * (1/p^2) ) = alpha/p + 1
(1/p -1) / ( v * (1/p^2) ) - 1 = alpha/p
p * ((1/p -1) / ( v * (1/p^2) ) - 1) = alpha

Let try:

alpha = 1

beta = 1

(p <- alpha/(alpha+beta))

(v <- (alpha*beta)/((alpha+beta)^2*(alpha+beta+1)))

(alpha_prime <- p * ((1/p -1 ) / ( v * (1/p^2) ) - 1) )

(beta_prime <-  alpha_prime * (1/p -1) )

It works