5  Probability Distributions - Beta

Open In Colab

Probability distributions describe the likelihood of different outcomes for a random event.

Here are a few ways to think about it:

5.1 Beta Distribution

Earlier we saw that binomial distribution is used to find the probability of successful outcomes, given a number of trials and the probability of the successful outcome.

In Beta distribution, we want to find the probability of success given we have observed the number of successful outcomes for a given number of trials.

For example we would use beta distribution to find the probability of flipping heads when we have observed 50 tosses of a coin and 30 of those were heads.

Code 
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import beta

# Parameters for the Beta distribution
alpha = 30
beta_param = 20 # Renamed to avoid conflict with the beta function

# Generate x-values (probabilities)
x = np.linspace(0, 1, 100) # Beta distribution is defined for x between 0 and 1

# Calculate the probability density function (PDF) for each x-value
pdf_values = beta.pdf(x, alpha, beta_param)

# Create the plot
plt.figure(figsize=(10, 6))
plt.plot(x, pdf_values)

# Add labels and title
plt.xlabel("Probability")
plt.ylabel("Density")
plt.title(f"Beta Distribution (alpha={alpha}, beta={beta_param})")
plt.grid(True)

# Display the plot
plt.show()

  • Beta distribution is different from binomial distribution as in binomial distribution we are looking at the distribution of k, the number of outcomes we care about, which is always something we can count.
  • For beta distribution, we are looking at the distribution of p, for which we have an infinite number of possible values.
  • The x axis of a beta distribution ranges from 0 to 1
  • The y axis is the probability density and can be thought of as the weight of each hypothesis.

It’s given by: \[Beta(p;α,β) = \frac{p^{α-1}(1-p)^{β-1}}{beta(α,β)}\]

  • p = probability of an event. This is our different hypothesis for the possible probabilities
  • α = how many times we observe an event we care about
  • β = how many times we observe the event we care about didn’t happen

The mean of a beta distribution is given by
\(μ = \frac{α}{(α + β)}\)
\(σ = \frac{(α*β)}{(α+β)^2*(α+β+1)}\)

If you have a sense of the average probability and the variance you can plug it into the equations above to get alpha and beta values.

5.1.1 Plotting Beta distributions for various alpha and beta values

If you know the mean and variance of a beta distribution you can also use bigger or smaller values of alpha and beta to get different distributions based on your understanding of data (ie prior beliefs). If you are uncertain you can choose small values, if you have stronger beliefs you can choose larger values and still keep the mean and variance same.

This flexibility allows you to incorporate your existing knowledge or lack thereof into the model.

Code 
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import beta

# Define a range of alpha and beta values to explore
alpha_values = [0.5,  2, 5, 10]
beta_values = [0.5, 2, 5, 10]

# Create a grid of subplots
fig, axes = plt.subplots(len(alpha_values), len(beta_values), figsize=(14, 10))

# Generate x-values (probabilities)
x = np.linspace(0, 1, 100)

# Iterate through the alpha and beta values and plot each distribution
for i, alpha in enumerate(alpha_values):
    for j, beta_param in enumerate(beta_values):
        # Calculate the probability density function (PDF)
        pdf_values = beta.pdf(x, alpha, beta_param)

        # Plot the distribution
        axes[i, j].plot(x, pdf_values)

        # Add titles and labels
        axes[i, j].set_title(f"α={alpha}, β={beta_param}")
        axes[i, j].set_xlabel("Probability")
        axes[i, j].set_ylabel("Density")
        axes[i, j].grid(True)

# Adjust layout to prevent overlapping titles and labels
plt.tight_layout()

# Display the plot
plt.show()

In real life situations, we almost never know the true probabilities for events.

That’s why the beta distribution is one of the most powerful tools for understanding data.

5.1.2 Is this coin fake?

Let’s say you are flipping a coin and notice that you get 12 heads and 18 tails and want to know the probability that the chance of getting heads is less than or equal to 0.5

import scipy.stats as stats
alpha = 12
beta = 18
threshold = 0.5
probability = stats.beta.cdf(threshold, alpha, beta)
print(probability)
0.867534551769495

OR you can also write this as

from scipy.integrate import quad
from scipy.stats import beta
# Define the function to integrate
def integrand(p):
 return beta.pdf(p, 12, 18)
# Perform the integration
result, error = quad(integrand, 0, 0.5)
result
0.8675345517694951

If you had to find the probability that the chance of getting heads was greater than 0.5 you would change the above code blocks like so

import scipy.stats as stats
alpha = 12
beta = 18
threshold = 0.5
probability = 1- stats.beta.cdf(threshold, alpha, beta)
print(probability)
0.132465448230505
from scipy.integrate import quad
from scipy.stats import beta
# Define the function to integrate
def integrand(p):
 return beta.pdf(p, 12, 18)
# Perform the integration
result, error = quad(integrand, 0.5, 1)
result
0.13246544823050507

5.1.3 Betting on Mbappe

Your friend has recently gotten into sports betting and loves watching Real Madrid. He wants to bet on Madrid forward Kylian Mbappe scoring a goal in the first game of the season and has asked for your help.

He is willing to bet only if there is at least a 0.5 probability that Mbappe’s scoring chance is greater than 70%. Should he make the bet?

You decide to pull up Mbappe’s stats and see that he has scored 332 goals in 462 games.
That makes α=332 and β=130

Code 
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import beta

# Parameters for the Beta distribution based on Mbappe's stats
alpha = 332
total_games = 462
beta_param = total_games - alpha

# Generate x-values (probabilities)
x = np.linspace(0, 1, 1000) # Beta distribution is defined for x between 0 and 1

# Calculate the probability density function (PDF) for each x-value
pdf_values = beta.pdf(x, alpha, beta_param)

# Create the plot
plt.figure(figsize=(10, 6))
plt.plot(x, pdf_values)

# Highlight the area under the curve from 0.7 to 1
x_fill = np.linspace(0.7, 1, 100)
pdf_fill = beta.pdf(x_fill, alpha, beta_param)
plt.fill_between(x_fill, pdf_fill, color='skyblue', alpha=0.5)


# Add labels and title
plt.xlabel("Probability")
plt.ylabel("Density")
plt.title(f"Beta Distribution for Mbappe's scoring chance (alpha={alpha}, beta={beta_param})")
plt.grid(True)

# Display the plot
plt.show()

Looking at the x axis you can see that all data is less than 0.8. We want to know how much exactly is above 0.7, the value our friend cares about for the bet (highlighted in blue above).

from scipy.stats import beta

# Parameters for the beta distribution based on Mbappe's stats
# Alpha is the number of successes (goals)
alpha = 332
# Beta is the number of failures (games without scoring)
total_games = 462
beta_param = total_games - alpha

# Calculate the integral from 0.7 to 1
# This is equivalent to 1 - CDF(0.7)
threshold = 0.7
probability = 1 - beta.cdf(threshold, a=alpha, b=beta_param)

print(f"The probability that Mbappe's chance of scoring is greater than {threshold*100}% is: {probability}")

# Check if this probability is greater than 0.5
bet_condition = 0.5
if probability > bet_condition:
    print(f"""Since the probability ({probability:.4f}) is greater than {bet_condition},
    you might consider making the bet.""")
else:
    print(f"""Since the probability ({probability:.4f}) is not greater than {bet_condition},
    he not should consider making the bet.""")
The probability that Mbappe's chance of scoring is greater than 70.0% is: 0.8140788834249968
Since the probability (0.8141) is greater than 0.5, 
    you might consider making the bet.