Продукт двух бета-распределений

Question

Скажем, у меня есть две случайные переменные:

X ~ Beta (α1, β1)

Y ~ Бета (α2, β2)

Я хотел бы вычислить распределение Z = XY (произведение случайных величин)

с scipy, я могу получить PDF единого Beta с:

from scipy.stats import beta 
rv = beta(a, b) 
x = np.linspace(start=0, stop=1, num=200) 
my_pdf = rv.pdf(x) 

Но как насчет продукта двух бета-версий? Могу ли я сделать это аналитически? (Решения Python/Julia/R в порядке).

Answer 1

аналитическое решение, посмотрите на this paper и this answer.

численный подход в R

set.seed(1) # for reproducability 

n <- 100000 # number of random variables 

# first beta distribution 
a1 <- 0.5 
b1 <- 0.9 
X <- rbeta(n, a1, b1) 

# second beta distribution 
a2 <- 0.9 
b2 <- 0.5 
Y <- rbeta(n, a2, b2) 

# calculate product 
Z <- X * Y 

# Have a look at the distributions 
plot(density(Z), col = "red", main = "Distributions") 
lines(density(X), lty = 2) 
lines(density(Y), lty = 2) 

Answer 2

FWIW, то же самое в Python

from scipy import stats 
import statsmodels.api as sm 
import matplotlib.pyplot as plt 

N = 100000 

y = stats.beta(.5, .9).rvs(N) 
x = stats.beta(.9, .5).rvs(N) 
z = x*y 
dens_z = sm.nonparametric.KDEUnivariate(z) 
dens_z.fit() 

dens_x = sm.nonparametric.KDEUnivariate(x) 
dens_x.fit() 

dens_y = sm.nonparametric.KDEUnivariate(y) 
dens_y.fit() 

fig, ax = plt.subplots() 
ax.plot(dens_z.support, dens_z.density, label='z') 
ax.plot(dens_x.support, dens_x.density, label='x') 
ax.plot(dens_y.support, dens_y.density, label='y') 
ax.legend() 
plt.draw_if_interactive()