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.

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

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

# 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")

# 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

# 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")

# 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{array}{r}\mathrm{Pr}[a\le X\le b]={\int }_a^b{f}_X(x)dx\end{array}$$

Properties:

1. Non-negativity: $f_X(x)\ge 0$ for all $x\in \mathbb{R}$.
2. Normalization: ${\int }_{-\mathrm{\infty }}^{\mathrm{\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{array}{r}{p}_X(k)=\mathrm{Pr}[X=k]\end{array}$$

Properties:

1. Non-negativity: $p_X(k)\ge 0$ for all integers $k$.
2. 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

# 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)" )

# 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

# 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" )

# 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{array}{r}{F}_X(x)=\mathrm{Pr}(X\le x)\end{array}$$

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{array}{r}{F}_X(x)={\int }_{-\mathrm{\infty }}^x{f}_X(u)du\end{array}$$

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 $-\mathrm{\infty }$ up to $x$.
• Since $\mathrm{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{array}{r}{f}_X(x)=\frac{d}{dx}{F}_X(x)\end{array}$$

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{array}{r}{F}_X(x)=\sum _{k\le x}{p}_X(k)\end{array}$$

• 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{array}{r}{p}_X(k)=F_X(k)-F_X(k^{-})\end{array}$$

  where $F_X(k^{-})=\underset{h\to 0^{+}}{lim}{F}_X(k-h)$ is the left-hand limit of $F_X$ at $k$.

Properties of the CDF

1. Non-decreasing: $F_X(x)$ is a non-decreasing function. If $x_1\le x_2$, then $F_X(x_1)\le F_X(x_2)$.

2. Limits:

   • Lower Bound:

     $$\begin{array}{r}\underset{x\to -\mathrm{\infty }}{lim}{F}_X(x)=0\end{array}$$

   • Upper Bound:

     $$\begin{array}{r}\underset{x\to \mathrm{\infty }}{lim}{F}_X(x)=1\end{array}$$

3. 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{array}{r}\underset{h\to 0^{+}}{lim}{F}_X(x+h)=F_X(x)\end{array}$$

   • 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{array}{r}\mathrm{Pr}(a<X\le b)=F_X(b)-F_X(a)\end{array}$$

  • Discrete Case:

    $$\begin{array}{r}\mathrm{Pr}(a<X\le b)=\sum _{k=a+1}^b{p}_X(k)\end{array}$$

• 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

# 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)" )

# 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

# 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" )

# 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{array}{r}{Q}_X(p)=F_X^{-1}(p)=inf\{x\in \mathbb{R}:F_X(x)\ge p\}\end{array}$$

• Discrete Case:

  $$\begin{array}{r}{Q}_X(p)=min\{k\in \mathbb{Z}:F_X(k)\ge p\}\end{array}$$

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

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

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

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

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{array}{r}{\mu }_k(c)=\mathrm{E}[(X-c{)}^k]\end{array}$$

• Raw Moments: When $c=0$, ${\mu }_k^{\prime }={\mu }_k(0)$.

• Central Moments: When $c=\mu$, the mean of $X$, ${\mu }_k={\mu }_k(\mu )$.

• Standardized Moments: ${\mu }_k^{\ast }={\displaystyle \frac{{\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

# 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

# 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

# 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

# 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

# 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

# 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

# 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

# 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

# Median in R median_continuous <- median(r_sample_continuous) median_continuous [1] 4.329341

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

Discrete Case

# Median in R median_discrete <- median(r_sample_discrete) median_discrete [1] 5

# 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.