Probability Functions in R and Python

Motivation

The intent of this post is to serve as a quick reference for working with basic probability functions in R and Python, covering both continuous and discrete cases. This guide provides some theoretical background, properties, and practical examples in both R and Python.

To demonstrate the equivalence between R functions and Python computations, we will use the same simulated samples in both languages throughout this guide. That is, the reticulate package in R will be used to pass data between R and Python.

Each code block has a slider to toggle between R and Python code snippets.

Packages and Libraries

We will utilize several packages in R and Python to perform our computations.

Show code

# Base R packages
library(stats)
library(moments)
library(EnvStats)
library(purrr)
library(dplyr)

Show code

from scipy.stats import gamma, chi2, binom, poisson, moment, describe
import numpy as np
import matplotlib.pyplot as plt

Simulate Samples

We begin by simulating samples from both continuous and discrete distributions, which will be used throughout this guide.

Continuous Distribution: Chi-squared Distribution

Show code

# Simulate samples in R from chi-squared distribution
r_sample_continuous <- rgamma(n = 1000, shape = 5 / 2, rate = 1 / 2)
hist(r_sample_continuous, main = "Histogram of Simulated Samples (R, Continuous)", xlab = "Value")

Show code

# Simulate samples in Python from chi-squared distribution
py_sample_continuous = gamma.rvs(a=5 / 2, scale=2, size=1000)
plt.hist(py_sample_continuous, bins=30)
plt.title("Histogram of Simulated Samples (Python, Continuous)")
plt.xlabel("Value")
plt.show()

Discrete Distribution: Binomial Distribution

Show code

# Simulate samples in R from binomial distribution
r_sample_discrete <- rbinom(n = 1000, size = 10, prob = 0.5)
hist(r_sample_discrete, breaks = seq(-0.5, 10.5, 1), main = "Histogram of Simulated Samples (R, Discrete)", xlab = "Value")

Show code

# Simulate samples in Python from binomial distribution
py_sample_discrete = binom.rvs(n=10, p=0.5, size=1000)
plt.hist(py_sample_discrete, bins=np.arange(-0.5,11.5,1))
plt.title("Histogram of Simulated Samples (Python, Discrete)")
plt.xlabel("Value")
plt.show()

Probability Density Function (Continuous) and Probability Mass Function (Discrete)

Probability Density Function (PDF) for Continuous Variables

The Probability Density Function (PDF) of a continuous random variable \(X\) is a function \(f_X(x)\) such that:

\[\begin{align*} \operatorname{Pr}[a \leq X \leq b] = \int_a^b f_X(x) \, dx \end{align*}\]

Properties:

Non-negativity: \(f_X(x) \geq 0\) for all \(x \in \mathbb{R}\).
Normalization: \(\int_{-\infty}^\infty f_X(x) \, dx = 1\).

Comments:

The PDF represents how densely the probability is distributed over the values of the random variable.
Although the probability at a single point is zero for the continous case, the PDF can be used to find the probability over an interval.

Probability Mass Function (PMF) for Discrete Variables

The Probability Mass Function (PMF) of a discrete random variable \(X\) is a function \(p_X(k)\) defined for integer values \(k\) such that:

\[\begin{align*} p_X(k) = \operatorname{Pr}[X = k] \end{align*}\]

Properties:

Non-negativity: \(p_X(k) \geq 0\) for all integers \(k\).
Normalization: \(\sum_{k} p_X(k) = 1\).

Comments:

The PMF gives the exact probability that the random variable takes on a specific value.
It is useful when the random variable can only take on discrete values.

Examples

Continuous Case: PDF of Chi-squared Distribution

Show code

# Compute PDF values in R
r_densities_continuous <- dchisq(r_sample_continuous, df = 10)
plot(r_sample_continuous, r_densities_continuous,
  xlab = "Quantiles", ylab = "Density",
  main = "PDF of Chi-squared Distribution (R)"
)

Show code

# Compute PDF values in Python
py_densities_continuous = chi2.pdf(r.r_sample_continuous, df=10)
plt.scatter(r.r_sample_continuous, py_densities_continuous);
plt.xlabel("Quantiles")
plt.ylabel("Density")
plt.title("PDF of Chi-squared Distribution (Python)")
plt.show()

Discrete Case: PMF of Binomial Distribution

Show code

# Compute PMF values in R
r_pmf_discrete <- dbinom(r_sample_discrete, size = 10, prob = 0.5)
plot(unique(r_sample_discrete), r_pmf_discrete[!duplicated(r_sample_discrete)],
  xlab = "Outcomes", ylab = "Probability",
  main = "PMF of Binomial Distribution (R)", type = "h"
)

Show code

# Compute PMF values in Python
py_pmf_discrete = binom.pmf(r.r_sample_discrete, n=10, p=0.5)
plt.stem(r.r_sample_discrete, py_pmf_discrete)
plt.xlabel("Outcomes")
plt.ylabel("Probability")
plt.title("PMF of Binomial Distribution (Python)")
plt.show()

Interpretation

Continuous Case: The PDF indicates the relative likelihood of the random variable near a point. Since the probability at an exact point is zero for continuous variables, we consider intervals.
Discrete Case: The PMF gives the exact probability of the random variable being equal to a specific value.

Cumulative Distribution Function (CDF)

The Cumulative Distribution Function (CDF) of a random variable \(X\) is a fundamental concept in probability theory that provides a complete description of the distribution of \(X\). It is defined as:

\[\begin{align*} F_X(x) = \operatorname{Pr}(X \leq x) \end{align*}\]

This definition holds for both continuous and discrete random variables. However, the way the CDF is calculated differs between the two cases due to the nature of their probability distributions.

Continuous Random Variables

For a continuous random variable, the CDF is calculated using an integral of its Probability Density Function (PDF) \(f_X(x)\):

\[\begin{align*} F_X(x) = \int_{-\infty}^{x} f_X(u) \, du \end{align*}\]

Key Points:

The PDF \(f_X(x)\) represents the density of probability at each point \(x\).
The CDF \(F_X(x)\) accumulates the probabilities from the lower bound \(-\infty\) up to \(x\).
Since \(\operatorname{Pr}(X = x) = 0\) for continuous variables, the CDF is a continuous and smooth function.

Relationship between PDF and CDF:

The PDF is the derivative of the CDF:

\[\begin{align*} f_X(x) = \frac{d}{dx} F_X(x) \end{align*}\]

Discrete Random Variables

For a discrete random variable, the CDF is calculated using a summation of its Probability Mass Function (PMF) \(p_X(k)\):

\[\begin{align*} F_X(x) = \sum_{k \leq x} p_X(k) \end{align*}\]

The PMF \(p_X(k)\) gives the probability that \(X\) takes the value \(k\).
The CDF \(F_X(x)\) sums up the probabilities for all possible values \(k\) that are less than or equal to \(x\).
The CDF is a step function, increasing at points where the random variable has positive probability mass.

Relationship between PMF and CDF:

The PMF can be obtained from the CDF by calculating the difference at each jump:

\[\begin{align*} p_X(k) = F_X(k) - F_X(k^-) \end{align*}\]

where \(F_X(k^-) = \lim_{h \to 0^+} F_X(k - h)\) is the left-hand limit of \(F_X\) at \(k\).

Properties of the CDF

Non-decreasing: \(F_X(x)\) is a non-decreasing function. If \(x_1 \leq x_2\), then \(F_X(x_1) \leq F_X(x_2)\).
Limits:
- Lower Bound:
  
  \[\begin{align*} \lim_{x \to -\infty} F_X(x) = 0 \end{align*}\]
- Upper Bound:
  
  \[\begin{align*} \lim_{x \to \infty} F_X(x) = 1 \end{align*}\]
Right-Continuity: The CDF is always right-continuous, regardless of whether the random variable is discrete, continuous, or mixed. This means that at any point \(x\):

\[\begin{align*} \lim_{h \to 0^+} F_X(x + h) = F_X(x) \end{align*}\]
- For Discrete Variables: The CDF has jumps at points where the random variable has positive probability mass, but it does not have discontinuities from the right.
- For Continuous Variables: The CDF is continuous everywhere, so it is trivially right-continuous.

Comments

Accumulated Probability: The CDF represents the total probability accumulated up to a certain value \(x\). It answers the question: “What is the probability that the random variable \(X\) is less than or equal to \(x\)?”
Calculating Probabilities Over Intervals:
- Continuous Case:
  
  \[\begin{align*} \operatorname{Pr}(a < X \leq b) = F_X(b) - F_X(a) \end{align*}\]
- Discrete Case:
  
  \[\begin{align*} \operatorname{Pr}(a < X \leq b) = \sum_{k = a+1}^{b} p_X(k) \end{align*}\]
Continuous Variables: The CDF is a smooth, continuous curve that increases from 0 to 1.
Discrete Variables: The CDF is a step function with jumps at the points where the random variable has non-zero probability.
Total Probability: The CDF approaches 1 as \(x\) approaches infinity, reflecting that the total probability over the entire space is 1.
Differentiation and Integration: For continuous random variables, the PDF and CDF are related through differentiation and integration, as previously noted.
Discontinuities: In the discrete case, the size of the jump at a point \(k\) in the CDF is equal to \(p_X(k)\), the probability mass at that point.

Examples

Continuous Case: CDF of Chi-squared Distribution

Show code

# Compute CDF values in R
r_cdf_continuous <- pchisq(r_sample_continuous, df = 10)
plot(r_sample_continuous, r_cdf_continuous,
  xlab = "Quantiles", ylab = "Cumulative Probability",
  main = "CDF of Chi-squared Distribution (R)"
)

Show code

# Compute CDF values in Python
py_cdf_continuous = chi2.cdf(r.r_sample_continuous, df=10)
plt.scatter(r.r_sample_continuous, py_cdf_continuous)
plt.xlabel("Quantiles")
plt.ylabel("Cumulative Probability");
plt.title("CDF of Chi-squared Distribution (Python)")
plt.show()

Discrete Case: CDF of Binomial Distribution

Show code

# Compute CDF values in R
r_cdf_discrete <- pbinom(r_sample_discrete, size = 10, prob = 0.5)
plot(r_sample_discrete, r_cdf_discrete,
  xlab = "Outcomes", ylab = "Cumulative Probability",
  main = "CDF of Binomial Distribution (R)", type = "s"
)

Show code

# Compute CDF values in Python
py_cdf_discrete = binom.cdf(r.r_sample_discrete, n=10, p=0.5)
plt.step(r.r_sample_discrete, py_cdf_discrete, where='post')
plt.xlabel("Outcomes")
plt.ylabel("Cumulative Probability")
plt.title("CDF of Binomial Distribution (Python)")
plt.show()

Comments

Continuous Case: The CDF smoothly increases from 0 to 1 as the quantile increases.
Discrete Case: The CDF increases in steps, reflecting the discrete nature of the variable.

Quantile Function

The Quantile Function is the inverse of the CDF and maps probabilities to quantiles.

Definitions

For a given probability \(p \in [0,1]\), the quantile function \(Q_X(p)\) is defined as:

Continuous Case:

\[\begin{align*} Q_X(p) = F_X^{-1}(p) = \inf\{ x \in \mathbb{R} : F_X(x) \geq p \} \end{align*}\]
Discrete Case:

\[\begin{align*} Q_X(p) = \min\{ k \in \mathbb{Z} : F_X(k) \geq p \} \end{align*}\]

Comments:

The quantile function tells us the value below which a certain percentage of data falls.
In the discrete case, it provides the smallest integer where the cumulative probability meets or exceeds the given probability.

Examples

Continuous Case

Show code

r_cdf_continuous <- pchisq(r_sample_continuous, df = 10)
# Quantile function in R
r_quantiles_continuous <- qchisq(p = r_cdf_continuous, df = 10)
# Verify
all.equal(r_quantiles_continuous, r_sample_continuous)

[1] TRUE

Show code

py_cdf_continuous = chi2.cdf(py_sample_continuous, df=10)
# Quantile function in Python
py_quantiles_continuous = chi2.ppf(py_cdf_continuous, df=10)
# Verify
np.allclose(py_quantiles_continuous, py_sample_continuous)

True

Discrete Case

Show code

r_cdf_discrete <- pbinom(r_sample_discrete, size = 10, prob = 0.5)
# Quantile function in R
r_quantiles_discrete <- qbinom(p = r_cdf_discrete, size = 10, prob = 0.5)
# Verify
all.equal(r_quantiles_discrete, r_sample_discrete)

[1] TRUE

Show code

py_cdf_discrete = binom.cdf(py_sample_discrete, n=10, p=0.5)
# Quantile function in Python
py_quantiles_discrete = binom.ppf(py_cdf_discrete, n=10, p=0.5)
# Verify
np.allclose(py_quantiles_discrete, py_sample_discrete)

True

Comments

The quantile function allows us to find thresholds corresponding to specific probabilities.
This is useful in statistical analyses, such as determining critical values.

Moments

Moments provide important characteristics of a probability distribution, such as its shape and spread.

Definitions

k-th Moment about a Point c:

\[\begin{align*} \mu_k(c) = \mathrm{E}\left[(X - c)^k\right] \end{align*}\]
Raw Moments: When \(c = 0\), \(\mu_k' = \mu_k(0)\).
Central Moments: When \(c = \mu\), the mean of \(X\), \(\mu_k = \mu_k(\mu)\).
Standardized Moments: \(\mu_k^* = \dfrac{\mu_k}{\sigma^k}\), where \(\sigma\) is the standard deviation.

Comments:

The first moment (mean) measures central tendency.
The second moment (variance) measures dispersion.
Higher moments (skewness and kurtosis) describe the shape of the distribution.

Computing Moments

Continuous Case

Show code

# First four central moments in R
central_moments_continuous_r <- map_dbl(1:4, ~ moment(r_sample_continuous, order = .x, central = TRUE))
names(central_moments_continuous_r) <- paste0("Moment_", 1:4)
central_moments_continuous_r

     Moment_1      Moment_2      Moment_3      Moment_4 
-4.733991e-16  1.080969e+01  4.478908e+01  6.202060e+02

Show code

# First four central moments in Python
central_moments_continuous_py = [moment(r.r_sample_continuous, moment=order) for order in range(1, 5)]
central_moments_continuous_py_dict = {f"Moment_{i+1}": m for i, m in enumerate(central_moments_continuous_py)}
central_moments_continuous_py_dict

{'Moment_1': np.float64(0.0), 'Moment_2': np.float64(10.80968784555429), 'Moment_3': np.float64(44.78908301891145), 'Moment_4': np.float64(620.2060122846481)}

Discrete Case

Show code

# First four central moments in R
central_moments_discrete_r <- map_dbl(1:4, ~ moment(r_sample_discrete, order = .x, central = TRUE))
names(central_moments_discrete_r) <- paste0("Moment_", 1:4)
central_moments_discrete_r

     Moment_1      Moment_2      Moment_3      Moment_4 
-5.675460e-16  2.447424e+00 -4.582284e-01  1.735153e+01

Show code

# First four central moments in Python
central_moments_discrete_py = [moment(r.r_sample_discrete, moment=order) for order in range(1, 5)]
central_moments_discrete_py_dict = {f"Moment_{i+1}": m for i, m in enumerate(central_moments_discrete_py)}
central_moments_discrete_py_dict

{'Moment_1': np.float64(0.0), 'Moment_2': np.float64(2.447424), 'Moment_3': np.float64(-0.4582283520000002), 'Moment_4': np.float64(17.351531292672)}

Comments

Moments help in understanding the distribution’s characteristics.
Comparing moments between continuous and discrete distributions can provide insights into their differences.

Descriptive Statistics

Descriptive statistics summarize key aspects of a dataset.

Computing Descriptive Statistics

Continuous Case

Show code

# Descriptive statistics in R
describe_continuous <- function(x) {
  list(
    nobs = length(x),
    minmax = range(x),
    mean = mean(x),
    var = var(x),
    skew = skewness(x),
    kurtosis = kurtosis(x, method = "fisher")
  )
}
describe_continuous(r_sample_continuous)

$nobs
[1] 1000

$minmax
[1]  0.1197989 24.5023743

$mean
[1] 5.073486

$var
[1] 10.82051

$skew
[1] 1.262132

$kurtosis
[1] 2.325345

Show code

# Descriptive statistics in Python
stats_continuous = describe(r.r_sample_continuous)
stats_continuous_dict = {
    'nobs': stats_continuous.nobs,
    'minmax': stats_continuous.minmax,
    'mean': stats_continuous.mean,
    'variance': stats_continuous.variance,
    'skewness': stats_continuous.skewness,
    'kurtosis': stats_continuous.kurtosis
}
for key, value in stats_continuous_dict.items():
    print(f"{key}: {value}")

nobs: 1000
minmax: (np.float64(0.11979889727592653), np.float64(24.502374324005203))
mean: 5.073485668442164
variance: 10.820508353908199
skewness: 1.2602376256130556
kurtosis: 2.307740310813246

Discrete Case

Show code

# Descriptive statistics in R
describe_discrete <- function(x) {
  list(
    nobs = length(x),
    minmax = range(x),
    mean = mean(x),
    var = var(x),
    skew = skewness(x),
    kurtosis = kurtosis(x, method = "fisher")
  )
}
describe_discrete(r_sample_discrete)

$nobs
[1] 1000

$minmax
[1]  0 10

$mean
[1] 5.024

$var
[1] 2.449874

$skew
[1] -0.1198589

$kurtosis
[1] -0.09768797

Show code

# Descriptive statistics in Python
stats_discrete = describe(r.r_sample_discrete)
stats_discrete_dict = {
    'nobs': stats_discrete.nobs,
    'minmax': stats_discrete.minmax,
    'mean': stats_discrete.mean,
    'variance': stats_discrete.variance,
    'skewness': stats_discrete.skewness,
    'kurtosis': stats_discrete.kurtosis
}
for key, value in stats_discrete_dict.items():
    print(f"{key}: {value}")

nobs: 1000
minmax: (np.int64(0), np.int64(10))
mean: 5.024
variance: 2.449873873873874
skewness: -0.11967905052475246
kurtosis: -0.10319421717102983

Computing the Median

Continuous Case

Show code

# Median in R
median_continuous <- median(r_sample_continuous)
median_continuous

[1] 4.329341

Show code

# Median in Python
median_continuous_py = np.median(r.r_sample_continuous)
median_continuous_py

np.float64(4.329341267889218)

Discrete Case

Show code

# Median in R
median_discrete <- median(r_sample_discrete)
median_discrete

[1] 5

Show code

# Median in Python
median_discrete_py = np.median(r.r_sample_discrete)
median_discrete_py

np.float64(5.0)

Comments

Descriptive statistics provide a summary of the dataset’s main characteristics.
They are crucial for initial data analysis and understanding the underlying distribution.