Density of the max of independent Beta random variables implementing https://math.stackexchange.com/a/2965587

demo.py

import pylab as plt
from math import prod
import numpy as np
from scipy.stats import beta as betarv
from scipy.special import betainc, beta as betafn

plt.ion()
plt.style.use('ggplot')
# %matplotlib qt

params = [(3, 4), (2, 3), (5, 2)]
betas = [betarv(a, b) for a, b in params]
size = 1_000_000

rvs = np.max(np.vstack([b.rvs(size) for b in betas]), axis=0)

xs = np.linspace(0, 1, 1001)


def pdf(x: float | int):
  return sum(b.pdf(x) * prod([betainc(*b2.args, x) for b2 in betas if b2 != b]) for b in betas)


plt.figure()
plt.hist(rvs, bins=50, density=True, label='histogram')
plt.plot(xs, [pdf(x) for x in xs], label='JimB PDF')
plt.xlabel('x')
plt.ylabel('P(x)')
plt.legend()
plt.savefig('hist-fit.png', dpi=300)
plt.savefig('hist-fit.svg')

hist-fit.svg