# Distribution of i.i.d. random variables $X$ and $Y$ if $XY \sim \text{Beta}(\alpha, \beta)$

This is a follow-up to the question Square root of a Beta(1,1) random variable, which received two great answers.

If $$XY \sim \text{Beta}(\alpha, \beta)$$, and $$X$$ and $$Y$$ are two independent identically distributed random variables, do we know their distribution? If there is no simple general solution, can we figure out the distribution of $$X$$ and $$Y$$ for the case $$XY \sim \text{Beta}(1, 1)$$?

From my search so far I guess this is a deconvolution problem for which no readily worked-out solution, such as the Kumaraswamy distribution in the previous question, exists. However, I would appreciate to be proven wrong in this regard.

In the case of $$Z=XY\sim \text{Beta}(\alpha,1)$$, the moment generating function (mgf) of $$-\ln(XY)=-\ln X-\ln Y$$ is \begin{align} M_{-\ln(XY)}(t) &= E(e^{-t\ln Z}) \\ &=E(Z^{-t}) \\ &=\int_0^1 z^{-t}\alpha z^{\alpha-1}dz \\ &= \frac{1}{1-t/\alpha} \end{align} which is the mgf of an exponentially distributed random variable with rate parameter $$\alpha$$.

Thus, if $$X$$ and $$Y$$ are iid, this means that $$-\ln X$$ and $$-\ln Y$$ must be Gamma distributed with shape parameter $$1/2$$ and rate parameter $$\alpha$$ (see wikipedia). Backtransforming, the pdf of $$X$$ and $$Y$$ is $$f(x)=\sqrt{-\frac\alpha{\pi\ln x}}x^{\alpha-1}$$ for $$0\le x\le 1$$.

A similar calculation in the general Beta$$(\alpha,\beta)$$-case leads to $$M_{-\ln(XY)}(t)=\frac{\Gamma(\alpha-t)\Gamma(\alpha+\beta)}{\Gamma(\alpha+\beta-t)\Gamma(\alpha)}$$ and the mgf of $$-\ln X$$ and $$-\ln Y$$ would need to be the square root of this. But that is not the mgf of any known distribution and such a distribution (that when convolved with itself produces the target density) may not even exist.

The pdf of $$-\ln X$$ can be computed via numerical integration of the inversion formula for the characteristic function $$\varphi(t)=(M_{-\ln(XY}(it))^{1/2}$$. This seems to produce valid pdfs for $$\beta\le 2$$ but for $$\beta > 2$$, the resulting function no longer appear to be a valid pdf: R code:

library(pracma)
f <- function(x,phi,...) {
integrand <- function(t) {
exp(-1i*t*x)*phi(t,...)
}
line_integral(integrand,c(-300+0i,300+0i))/(2*pi)
}
f <- Vectorize(f,vectorize.args = "x")
phi <- function(t, alpha, beta) {
sqrt(gammaz(alpha - 1i*t)*gammaz(alpha + beta)/(gammaz(alpha + beta - 1i*t)*gamma(alpha)))
}
x <- seq(-1,3,by=.01)
fx <- f(x, phi, alpha=3,beta=2.5)
plot(x,Re(fx),type="l")

• This looks extremely useful. Can you maybe add a hint, how you derived the relationships in the first paragraph? Or how they are called? May 4 at 9:08
• @LuckyPal I added some more explanation and made the result slightly more general. May 4 at 9:37
• "but for β>2, the resulting function no longer appear to be a valid pdf" Is it not valid because of the negative values? But, those negative values only appear outside of the support $0<x<1$. 7 hours ago