Calculate percentiles of probabilistic functions

Question

For easier understanding, let's assume we talk about sport betting. There are 3 independent games with team_A winning with a probability of 70%, 50%, 10% respectively. Each bet has its odds (odds_1, odds_2, odds_3). If was to bet a certain amount in each of the games (a_1, a_2, a_3), here is my gain/loss:

(a_1 * odds_1 with 70% prob. else -a_1) + 
(a_2 * odds_2 with 50% prob. else -a_2) + 
(a_3 * odds_3 with 10% prob. else -a_3)

In the worst case, I will lose (a_1 + a_2 + a_3). In the best case, I will win (a_1*odds_1 + a_2*odds_2 + a_3*odds_3).

Is there a way to calculate the percentiles of the expected gains/losses assuming the probabilities are the actual ones? I could do that by running many simulations and get the estimated percentiles. I'm wondering if there is an analytical way to estimate the percentiles without simulation.

As a concrete example, odds could be (1.3, 1.9, 4.3), the probabilities (70%, 50%, 10%) and the amounts (0.5, 1.2, 0.25).

Edit: In practice, I have from 60 to 80 such cases so there are many possible outcomes. Moreover, I want to be able to calculate these percentiles for thousands of different amounts combinations. Therefore, exhaustive calculation is not possible.

Answer 1

As suggested in an answer to a similar question, a brute-force calculation can do fine.

There's a trick, though: when the rewards $a_i$ and probabilities $p_i$ are arbitrary, each bet has two distinct possible gains ($a_i(1-p_i)/p_i$ and $b_i = -a_i$) which implies the total gain after $n$ bets could be a number in a set of up to $2^n$ distinct values, which is huge. So round them. Although the rounding creates a discretization error, in the sums those errors will tend to wash out and you can control the maximum error by discretizing the problem more finely.

Because we're rounding and the problem doesn't change when we change the units of currency, we may as well assume all rewards have been rounded to the nearest whole number. A probability distribution of rewards can therefore be expressed as a vector

$$\mathbf{q} = \ldots q_{-2},q_{-1}, q_0, q_1, q_2, q_3, \ldots$$

where $q_n$ is the chance that the sum of your rounded gains (and losses) equals $n.$ For indices $n$ less than the sum of all possible losses $-B$ or greater than the sum of all possible gains $A$, $q_n$ must be zero, so this is really a finite vector.

When you bet again, with possible gain $a$ (with probability $p$) or loss $-b$ (with probability $1-p$), the axioms of probability imply $\mathbf q$ is updated to

$$\mathbf{q}^\prime = (1-p) \mathcal{S}_{-b}[\mathbf{q}] + p\, \mathcal{S}_a[\mathbf{q}]$$

where $\mathcal{S}_i$ shifts a vector by $i$ places (to the left for negative $i$ and to the right for positive $i$):

$$(\mathcal{S}_i[\mathbf q])_j = q_{j-i}$$

for all $j.$ Because $\mathbf q$ has nonzero entries confined to $A+B+1$ consecutive locations, computing this update takes $O(A+B)$ computational effort. Starting with $\mathbf{q} = (\ldots, 0, 0, 1, 0, 0, \ldots)$ (with $q_0=1$), $n$ of these updates will produce the answer with $O(n(A+B))$ effort. Because a single core can perform hundreds of millions of basic operations a second in R (and even more in C++ or Python, for instance), the complete distribution--not just its quantiles--can be made available in a fraction of a second.

This can be appreciably sped up by periodically dropping tiny probabilities at the tails of the distribution.

In the crude R implementation shown below, which has no optimization, a problem with $n=8000$ bets discretized into approximately $80,000$ bins was computed in four seconds. That the resulting distribution used only $1188$ bins suggests the preceding optimization could speed it up by orders of magnitude--but it all depends on the patterns of bets and their odds. When there are a tiny number of bets with huge risks (very large $a$ and $b,$) the range of total rewards that have an appreciable chance of occurring can be large.

To illustrate, I randomly generated $80$ bets having widely varied odds and payoffs. As described in the question, all bets were constructed to have expected values of zero. Here is a scatterplot matrix (on log-log scales) of the data:

You can see the wide variation in payoffs and probabilities.

Next is a plot of the distribution function of the payoffs. Remember, because rewards were rounded, this is a discrete distribution: it's really a bar plot with so many thin bars they crowd into a shaded region. The bin widths are $1.08:$

Total computation time was 0.05 seconds.

This is a problematic distribution in that it has a long right tail: this reflects some small prospects of making huge gains. It also demonstrates that the Central Limit Theorem does not necessarily provide any insight into the solution, as some people might initially hope.

Finally, I checked this calculation (and many others like it) by simulating the sum of these 80 bets one million times. (That took 17 seconds, some 300 times longer than the calculation.) I computed 1st, 2nd, ..., 99th percentiles of the simulated data and graphed them on this QQ plot against the corresponding percentiles of the computed distribution:

The QQ plot is drawn in black and over it is plotted in red the line of equality: they match almost exactly at this resolution, demonstrating close agreement.

---

Given the probability vector p, positive payoff vector a, and loss vector b, this code ultimately computes the quantile function f, which can be applied to any value between $0$ and $1$ to return (very quickly) the corresponding quantile of the distribution of total payoffs from the bets.

#
# The calculation.
# The entire distribution is returned in a list.
# The result will be discretized into about 2*N bins.
#
binary.sum <- function(p, a, b, N=1000, thresh=1e-8) {
  i <- order(a+b) # (It's more efficient to process the smaller bets first)
  p <- p[i]
  d <- (b+a)[i]

  x.max <- sum(a)
  x.min <- -sum(b)
  h <- (x.max - x.min) / (2*N)   # Common bin width
  ix <- function(x) round(x / h) # Assign a value to a bin

  M <- sum(ix(d)) + 1
  z <- c(1, rep(0, M-1))
  n <- 1 # Logical length of z
  for (i in seq_along(p)) {
    y <- z[1:n]
    k <- ix(d[i])
    z[1:n] <- (1-p[i]) * y
    z[1:n + k] <- p[i] * y + z[1:n + k]
    n <- n + k
  }
  x <- seq(x.min, M*h+x.min, length.out=M+1)
  i.pos <- z > thresh # Zero out any probabilities of `thresh` or less
  i <- which.max(i.pos)
  j <- M+1 - which.max(rev(i.pos))
  list(start=x[i], stop=x[j], h=h, prob=z[i:j], x=x[i:j])
}
# Create an example.
#
set.seed(17)
n <- 80
p <- (1/n + rgamma(n, 1/4))
p <- p / (1+max(p))                    # The probabilities
b <- rgamma(n, 1, 1/20)                # A "success" nets `a`; a failure loses `b`.
a <- b * (1-p) / p
#-- Figure 1:
pairs(cbind(Probability=p, `Gain (a)`=a, `Loss (b)`=b), log="xy")
system.time(
  q <- binary.sum(p, a, b, N=1e5, thresh=1e-6)
)
#
# Figure 1.
#
with(q, plot(x, prob, type="h", xlab="Total Payoff", ylab="Probability",
     main="Computed Distribution"))
#
# Simulation.
#
system.time({
  sim <- replicate(1e6, sum(ifelse(rbinom(length(p), 1, p)==1, a, -b)))
})
#-- Optional ECDF plot of the simulated results on which the computed CDF is superimposed.
if (length(sim) <= 1e4) {
  plot(ecdf(sim), verticals=TRUE, lwd=2, main="Simulated ECDF with Computed CDF")
  lines(x, cumsum(q$prob), type="s", lwd=2, col="Red")
}
#
# Check the quantile function with a QQ plot.
#
f <- with(q, stepfun(cumsum(prob), c(start-h, x))) # The computed quantile function
qq <- seq(0, 1, length.out=101)
qq <- qq[2:(length(qq)-1)]
#-- Figure 3:
plot(f(qq), type="l", quantile(sim, qq), asp=1, lwd=3,
     main="QQ Plot",
     xlab="Computed quantile", ylab="Quantile in simulation")
abline(0:1, col="#e81010c0", lwd=2)

Answer 2

Let $p_i$ be the probability of winning game $i$, and $X_i\sim\operatorname{Bernoulli}(p_i)$ be the indicator variable of winning game $i$. If we bet $a_i$ dollars on game $i$ with odds $o_i$, the final amount of money we have is $$\sum_{i=1}^n X_i(a_io_i + a_i) - \sum_{i=1}^n a_i.$$ Let $Y_i = X_i(a_io_i + a_i) - a_i$, so $Y_i$ has mean $\mu_i = p_ia_io_i + p_ia_i - a_i$ and variance $\sigma_i^2 = (a_io_i + a_i)^2p(1-p)$. For small $n$, we can directly compute the probabilities and total earnings for each of the $2^n$ possible outcomes. For large $n$, we can use Lindeberg's CLT which tells us that $$\frac{\sum_{i=1}^n (Y_i - \mu_i)}{\sqrt{\sum_{i=1}^n \sigma_i^2}}\approx N(0, 1).$$ By this approximation, we have $$\sum_{i=1}^n Y_i\approx N\left(\sum_{i=1}^n (p_i a_io_i + p_ia_i - a_i), \sum_{i=1}^npa_io_i + pa_i - a_i\right).$$

Of course, how well this works on your data, you'll have to test that directly. I would recommend running some simulations to see how well this theoretical approximation estimates the true data (as it may be the case that you don't have a large enough $n$ for which CLT may be accurate).