Expected value of the max of scaled powers of the same random variable

Question

I've been looking at order statistics and the behavior of expectations when max is involved, but that literature always discussed iid random variables whereas I have a strange situation where I have only one random variable that's being scaled and exponentiated by deterministic quantities so it's not clear to me if that's relevant.

I have a random variable $X \sim Beta(\alpha, \beta)$, and I want to compute

$$\mu = E[\max \lbrace a_1 X^{b_1}, a_2 X^{b_2}, \cdots, a_n X^{b_n} \rbrace]$$

for known $a_i$ and $b_i$. If it helps, my $a_i \in (0, 1]$ are strictly increasing and $b_i \in (0, 1]$ also strictly increasing.

I know that $E[X^c] = \frac{B(\alpha + c, \beta)}{B(\alpha, \beta)}$ (where $B(\cdot, \cdot)$ is the Beta function), and I can convince myself numerically that if all $a_i = 1$, then $\mu = \max_i \lbrace E[X^{b_i}] \rbrace$. But I'm not sure how to show this analytically so I'm not sure what to do when I see that numerically, $\mu \neq \max_i \lbrace a_i E[X^{b_i}] \rbrace$, there's a small but regular bias in my Monte Carlo simulation.

I feel like this case, where I just have a single random variable, should be quite straightforward to handle but obviously I'm missing some perspective in handling situations like this. Can I make progress with this?

---

Some Python code leveraging Numpy and Scipy to reproduce the results above with some example numbers follows—

from scipy.stats import beta as betarv
from scipy.special import gamma
import numpy as np


def monteCarloTest(alpha=2, beta=7, size=1_000_000, n=5):
  xs = betarv.rvs(alpha, beta, size=size)

  bvec = (np.arange(n + 1) / (n + 1 - 1))[1:]
  arr = xs**(bvec[:, np.newaxis])
  monteCarloA1 = np.mean(np.max(arr, axis=0))

  analyticalMeans = gamma(alpha + beta) * gamma(alpha + bvec) / (
      gamma(alpha) * gamma(alpha + beta + bvec))
  analyticalA1 = np.max(analyticalMeans)
  # Ok that makes sense, E[max(X ** b)] = max(E[X ** b]), here b is known deterministic vector and X is Beta rv

  # what about E[max(a * (X ** b))] for a and b known vectors?
  avec = np.linspace(.1, .9, n)
  monteCarlo = np.mean(np.max(avec[:, np.newaxis] * arr, axis=0))
  analytical = np.max(analyticalMeans * avec)
  print(f'{monteCarloA1=:.5f}, {analyticalA1=:.5f}; {monteCarlo=:.5f}, {analytical=:.5f}')


monteCarloTest()
monteCarloTest()
monteCarloTest()
monteCarloTest()
monteCarloTest()
monteCarloTest()

Running the above produces the following output for me:

monteCarloA1=0.71648, analyticalA1=0.71640; monteCarlo=0.21708, analytical=0.20414
monteCarloA1=0.71641, analyticalA1=0.71640; monteCarlo=0.21696, analytical=0.20414
monteCarloA1=0.71653, analyticalA1=0.71640; monteCarlo=0.21703, analytical=0.20414
monteCarloA1=0.71646, analyticalA1=0.71640; monteCarlo=0.21705, analytical=0.20414
monteCarloA1=0.71639, analyticalA1=0.71640; monteCarlo=0.21689, analytical=0.20414
monteCarloA1=0.71637, analyticalA1=0.71640; monteCarlo=0.21690, analytical=0.20414

The first two numbers on each line are the Monte Carlo vs "analytical" expectation for all $a_i=1$; these are in tight agreement which is what leads me to my belief that the max can be moved into the expectation for this simpler case.

The second two numbers have $a_i$ varying, and give the Monte Carlo and my straightforward-but-incorrect analytical expectation. I note the latter is persistently below the former, indicating the mismatch.

Answer 1

I just needed to look at a plot of the $n$ functions going into $\max$ to get a sense of how to approach this:

Above, I added a dotted black line representing the max, and you see it follow the first $i=1$ blue line, then switch to the second $i=2$ orange line, and then the $i=3$ green line (and so on to the rest). So: $\max$ is choosing one of the $n$ functions over a specific domain and I know the $n-1$ crossing points between the $i$th function $a_i x^{b_i}$ and the $i+1$th: some algebra to solve for $a x^{b} = c x^{d}$ leads to this expression for the crossings:

$$x_{i,i+1} = \left(\frac{a_{i+1}}{a_i}\right)^{\frac{1}{b_i - b_{i+1}}}$$

So the expectation is broken into $n$ integrals over contiguous non-overlapping intervals, $E\left[\max_i \left\lbrace a_i X^{b_i} \right\rbrace\right]$ is

$$\frac{1}{B(\alpha,\beta)}\int_0^{x_{1,2}}a_1 x^{b_1} x^{\alpha - 1}(1-x)^{\beta-1} dx + \\ \frac{1}{B(\alpha,\beta)}\int_{x_{1,2}}^{x_{2,3}}a_2 x^{b_2} x^{\alpha - 1}(1-x)^{\beta-1} dx + \\ \vdots \\ \frac{1}{B(\alpha,\beta)}\int_{x_{n-1,n}}^{1}a_n x^{b_n} x^{\alpha - 1}(1-x)^{\beta-1} dx$$

Given that

$$\frac{a}{B(\alpha,\beta)}\int_{x_1}^{x_2} x^{\alpha + b - 1} (1-x)^{\beta-1} dx = \frac{a B(\alpha+b, \beta)}{B(\alpha, \beta)} (B_{x_2}(\alpha+b,\beta) - B_{x_1}(\alpha+b,\beta))$$

where $B_c(\cdot,\cdot)$ is the regularized incomplete Beta function, I could adapt the Python code in my original question to the following to conveniently compute the expectation analytically. Now it matches the Monte Carlo simulation quite well!

---

from scipy.stats import beta as betarv
from scipy.special import gamma, betainc, beta as betafn
import numpy as np
from itertools import pairwise
from functools import cache


def monteCarloTest(alpha=2, beta=7, size=1_000_000, n=5, viz=False):
  xs = betarv.rvs(alpha, beta, size=size)

  bvec = (np.arange(n + 1) / (n + 1 - 1))[1:]
  arr = xs**(bvec[:, np.newaxis])
  monteCarloA1 = np.mean(np.max(arr, axis=0))

  analyticalMeans = gamma(alpha + beta) * gamma(alpha + bvec) / (
      gamma(alpha) * gamma(alpha + beta + bvec))
  analyticalA1 = np.max(analyticalMeans)
  # Ok that makes sense, E[max(X ** b)] = max(E[X ** b]), here b is known deterministic vector and X is Beta rv

  # what about E[max(a * (X ** b))] for a and b known vectors?
  avec = np.linspace(.1, .9, n)
  arr = avec[:, np.newaxis] * arr
  #   print(np.mean(arr, axis=1))
  #   print(avec * analyticalMeans)

  monteCarlo = np.mean(np.max(arr, axis=0))
  analytical = np.max(analyticalMeans * avec)

  # above: same code as question

  # below: new code: plotting and correct analytical expectation
  if viz:
    import pylab as plt
    plt.ion()
    xvec = np.linspace(0, 1, 10001)
    arr = avec[:, np.newaxis] * xvec**(bvec[:, np.newaxis])
    plt.figure()
    plt.semilogy(xvec, arr.T)
    plt.semilogy(xvec, np.max(arr, axis=0), 'k:')
    plt.xlabel('x')
    plt.ylabel('a_i x^b_i')

  crossings = [0.0] + [
      (a2 / a1)**(1 / (b1 - b2)) for ((a1, b1), (a2, b2)) in pairwise(zip(avec, bvec))
  ] + [1.0]

  betaincLoHi = lambda a, b, c1, c2: (cachedBetainc(a, b, c2) - cachedBetainc(a, b, c1)
                                     ) * cachedBeta(a, b)
  bab = cachedBeta(alpha, beta)
  correct = sum(a * betaincLoHi(alpha + b, beta, lo, hi)
                for (a, b, (lo, hi)) in zip(avec, bvec, pairwise(crossings))) / bab
  print(f'{monteCarloA1=:.5f}, {analyticalA1=:.5f}; {monteCarlo=:.5f}, {correct=:.5f}')


@cache
def cachedBetainc(*args, **kwargs):
  return betainc(*args, **kwargs)


@cache
def cachedBeta(*args, **kwargs):
  return betafn(*args, **kwargs)


monteCarloTest()
monteCarloTest()
monteCarloTest()
monteCarloTest()
monteCarloTest()
monteCarloTest()

(Caching the output of betafn and betainc is a helpful but naive optimization, makes this 5x faster on my benchmark, but obviously more analytical massaging will be better.)

This now produces the following output:

monteCarloA1=0.71657, analyticalA1=0.71640; monteCarlo=0.21711, correct=0.21695
monteCarloA1=0.71639, analyticalA1=0.71640; monteCarlo=0.21696, correct=0.21695
monteCarloA1=0.71613, analyticalA1=0.71640; monteCarlo=0.21665, correct=0.21695
monteCarloA1=0.71644, analyticalA1=0.71640; monteCarlo=0.21701, correct=0.21695
monteCarloA1=0.71643, analyticalA1=0.71640; monteCarlo=0.21697, correct=0.21695
monteCarloA1=0.71633, analyticalA1=0.71640; monteCarlo=0.21688, correct=0.21695