Packages and Library

# Should come with Base R
library(stats)
library(moments)
library(EnvStats)
library(purrr)
library(dplyr)

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

Simulate Samples

r_sample <- rgamma(n = 1000, shape = 5 / 2, rate = 1 / 2) hist(r_sample)

# The rate parameter in R is the inverse of the scale parameter 1 / scale in Python py_sample = gamma.rvs(a=5 / 2, scale=2, size=1000) plt.hist(py_sample); plt.show()

Probability Density Function

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

where \(f_X(x): \mathbb{R} \rightarrow[0,1]\).

Python

In Scipy, the pdf method provides the densities for the given array of quantiles.

r_densities <- dchisq(r_sample, df = 10) plot(r_sample, r_densities, xlab = "quantiles")

py_densities = chi2.pdf(r.r_sample, df=10) plt.scatter(r.r_sample, py_densities); plt.show()

---

Cumulative Distribution Function

\[\begin{align*} F_X(x)=\operatorname{Pr}(X \leq x)=\int_{-\infty}^x f_X(u) d u \end{align*}\]

where \(F_X(x): \mathbb{R} \rightarrow[0,1]\).

Python

In Scipy, the cdf method provides the probabilities for the given array of quantiles.

r_probabilities <- pchisq(q = r_sample, df = 10) plot(r_sample, r_probabilities)

py_probabilities = chi2.cdf(r.r_sample, df=10) plt.scatter(r.r_sample, py_probabilities); plt.show()

---

Quantile function

Let \(X\) be a random variable with the cumulative distribution function (CDF) \(F_X(x)\). Then, the function \(Q_X(p):[0,1] \rightarrow \mathbb{R}\), which is the inverse \(\mathrm{CDF}\), is the quantile function of \(X\). More precisely, the QF is the function that, for a given probability \(p \in[0,1]\), returns the smallest \(x\) for which \(F_X(x)=p\)

\[\begin{align*} Q_X(p)=\min \left\{x \in \mathbb{R} \mid F_X(x)=p\right\} \end{align*}\]

r_quantiles <- qchisq(p = r_probabilities, df = 10)
# The returned quantities should be equivalent to the simulated sample
all(near(r_quantiles, r_sample))

## [1] TRUE

py_quantiles = chi2.ppf(py_probabilities, df=10)
# The returned quantities should be equivalent to the simulated sample
np.allclose(py_quantiles, r.r_sample)

## True

Moments

Let \(X\) be a random variable, let \(c\) be a constant and let \(k\) be a positive integer. Then, the \(k\) -th moment of \(X\) about \(c\) is defined as the expected value of \(X\) minus \(c\) raised to the \(k\) -th power:

\[\begin{align*} \mu_k(c)=\mathrm{E}\left[(X-c)^k\right] \end{align*}\]

The \(k\)-th moment of \(X^k\) can be expressed as:

• the \(k\)-th raw moment \(\mu_k^{\prime}=\mu_k(0)\);
• the \(k\)-th central moment \(\mu_k=\mu_k(\mu)\);
• the \(k\)-th standardized moment \(\mu_k^*=\mu_k / \sigma^k\).

The k-th central moment of a data sample is:

\[\begin{align*} \mu_k(\mu)=\frac{1}{n - 1} \sum_{i=1}^n\left(x_i-\bar{x}\right)^k \end{align*}\]

where

• n is the number of samples
• \(\bar{x}\) is the sample mean

Note: we use \(n - 1\) since we lose one degree of freedom from using the sample mean to estimate the population mean.

# First four central moments
map2_dbl(
  .x = c(1, 2, 3, 4),
  .y = c(T, T, T, T),
  .f = ~ moment(x = r_sample, order = .x, central = .y)
)

## [1] -4.174439e-16  9.868398e+00  3.440686e+01  4.482257e+02

# First four central moments
# This implementation does not correct for degree of freedom
[moment(a=r.r_sample, moment=order) for order in [1, 2, 3, 4]]

## [np.float64(0.0), np.float64(9.868397946961613), np.float64(34.40685586626528), np.float64(448.2256977917112)]

Desciptive Statistics

• Sample Mean: First raw moment.

• Sample Median: If \(n=2m+1\) for some integer \(m\), the sample median is \(X_{(m+1)}\) and an order statistic. When \(n\) is even, \(n=2m\), there are two middle values, \(X_{(m)}\) and \(X_{(m+1)}\), and the sample median is some function of the two (usually the average) and hence not an order statistic.

• Variance: Second central moment.

• Skewness: Third standardized moment.

• Kurtosis: Fourth standardized moment.

describe <- function(x) {
  list(
    nobs = length(x),
    minmax = c(min(x), max(x)),
    mean = mean(x),
    var = var(x),
    skew = skewness(x),
    kurtosis = kurtosis(x, method = "fisher")
  )
}
describe(r_sample)

## $nobs
## [1] 1000
## 
## $minmax
## [1]  0.1183308 21.0500038
## 
## $mean
## [1] 5.08348
## 
## $var
## [1] 9.878276
## 
## $skew
## [1] 1.111545
## 
## $kurtosis
## [1] 1.616668

# Instance of 'scipy.stats.stats.DescribeResult'
results = {key:value for key, value in zip(describe(r.r_sample)._fields, describe(r.r_sample))}
for key in results:
    print(key, ' : ', results[key])

## nobs  :  1000
## minmax  :  (np.float64(0.11833077848532816), np.float64(21.050003760803623))
## mean  :  5.083480371960868
## variance  :  9.878276223184798
## skewness  :  1.109877396544084
## kurtosis  :  1.6026022291528061

median(r_sample)

## [1] 4.52105

np.median(r.r_sample)

## np.float64(4.5210499170500515)