Моделирование распределений на языке Python

Моделирование будем делать с помощью библиотеки scipy.

Полный список доступных распределений можно посмотреть по здесь: https://docs.scipy.org/doc/scipy/reference/stats.html

Далее на нормальном распределении будут приведены примеры работы с библиотекой.

1) Промоделируем случайную величину
2) Построим гистограмму и плостность распределения
3) Нарисуем графики функции распределения и накопленных частот

Для начала импортируем библиотеку:

In [1]:
from scipy import stats

Непрерывные распределения

Для моделирования распредения воспользуемся функцией .rvs()
Промоделируем одномерное нормальное распределение $N(0, 1)$:

In [2]:
normal = stats.norm(loc=0, scale=1)
sample = normal.rvs(size=1000, random_state=42)

Теперь, в массиве sample у нас находится выборка из нормального распределения, размера $1000$.

Построим гистограмму, для этого воспользуемся библиотекой matplotlib:

In [3]:
import matplotlib.pyplot as plt
%matplotlib inline
In [4]:
plt.hist(sample, bins=20) # Параметр bins - отвечает за количество разбиений в гистограмме
plt.title('Гистограмма выборки из нормального распределения')
plt.show()

Поверх гистограммы построим график функции плотности распределения, чтобы получить значение плотности распределения в заданной точке (или наборе точек), воспользуемся функцией pdf():

In [5]:
import numpy as np

x = np.linspace(-5, 5, 100) # Отрезок на котором будем строить график
plt.hist(sample, density=True, bins=20, label='histogram')
plt.plot(x, normal.pdf(x), 'r-', lw=3, label='frozen pdf')
plt.legend()
plt.show()

Построим графики накопленных частот и функции распределения, для этого воспользуемся функцией ECDF библиотеки statsmodels:

In [6]:
from statsmodels.distributions.empirical_distribution import ECDF

ecdf = ECDF(sample)

x = np.linspace(min(sample), max(sample))
plt.step(x, ecdf(x), label='Empirical CDF')
plt.plot(x, normal.cdf(x), label='Theoretical CDF')
plt.legend()
plt.show()

Давайте убедимся, что гистограмма при большом объеме выборки почти совпадает с теоретической плотностью

Для этого сгенерируем выборку размером 10000 из стандартного нормального распределения, чтобы это было наглядно, увеличим количество разбиений в гистограмме:

In [7]:
sample = stats.norm(loc=0, scale=1).rvs(size=10000, random_state=42)
x = np.linspace(-5, 5, 100) # Отрезок на котором будем строить график
plt.hist(sample, density=True, bins=50, label='histogram')
plt.plot(x, normal.pdf(x), 'r-', lw=3, label='frozen pdf')
plt.legend()
plt.show()

Дискретные распределения

Для моделирования распредения воспользуемся функцией .rvs()
Промоделируем распределение Пуассона с параметром $\lambda = 1$:

In [8]:
lmbda = 1
pois = stats.poisson(lmbda)
sample = pois.rvs(size=1000, random_state=42)

Теперь, в массиве sample у нас находится выборка из распределения Пуассона, размера $1000$.

Построим гистограмму:

In [9]:
values, counts = np.unique(sample, return_counts=True)

В массиве values находятся отсортированные уникальные значения смодедированной выборки, а в массиве counts -- количество повторений элементов из массива values.

In [10]:
plt.vlines(values, 0, counts, color='r', lw=4)
plt.title('Emperical Poisson Probability Mass Function (PMF)')
plt.show()

Добавим на график теоретические значения функции вероятности:

In [11]:
plt.vlines(values, 0, counts/sum(counts), lw=3, color='r', label='Emperical Poisson PMF')
plt.vlines(values, 0, pois.pmf(values), colors='k', linestyles='dashdot', lw=2, label='Poisson PMF')
plt.legend()
plt.show()
In [12]:
ecdf = ECDF(sample)

x = np.linspace(min(sample), max(sample))
plt.plot(x, ecdf(x), color='r', label='Empirical CDF')
plt.plot(x, pois.cdf(x), label='Theoretical CDF')
plt.legend(loc='lower right')
plt.show()

Стандартизация

Иногда требуется стандартизировать дискретное распределение (например для того, чтобы более осмысленно сравнивать разные по природе выборки или, например, "похожесть" дискретного распределения на нормальное). Это делается несложным преобразованием: вычитанием среднего и делением на стандартное отклонение

$$\frac{X_i - \bar x}{s},$$

где $X_i$ - элемент из выборки $\{X_i\}_{i = 1}^n,$ $\bar x$ - среднее, а $s$ - стандартное отклонение.

Промоделируем выборку из Пуассоновского распределения с достаточно большим параметром $\lambda = 100.$

In [13]:
lmbda = 100
pois = stats.poisson(lmbda)
sample = pois.rvs(size=100000, random_state=42)
_, bins, _nvm = plt.hist(sample, density=True, bins=15, label='histogram') 
# в переменной bins теперь разбиения гистограммы

plt.show()

Допустим, мы хотим сравнить это дискретное распределение со стандартным нормальным распределением, хотя бы визуально. Для этого нужно стандартизировать нашу выборку и построить гистограмму вместе с плотностью нормального распределения.

In [14]:
stand_sample = (sample - sample.mean())/ sample.std()

Заметин, что вместе со стандартизацией выборки, мы также применим тоже преобразование и к вектору bins, чтобы сохранить ширину и позицию столбцов гистограммы при стандартизации

In [15]:
new_bins = (bins - sample.mean()) / sample.std()
In [16]:
x = np.linspace(-4, 4, 100)
plt.hist(stand_sample, density=True, bins=new_bins)
plt.title('Гистограмма стандартизированной выборки и плотность N(0, 1)')
plt.plot(x, stats.norm.pdf(x), 'r-', lw=3)
plt.show()

В данном случае можно сделать интересное наблюдение, которое соответсвует действительности: при больших $\lambda$ Пуассоновское распределение неплохо аппроксимируется нормальным распределением со средним и диспресией, равными $\lambda$