Физтех.Статистика
Скачать ipynb
Python для анализа данных¶
Работа со случайными величинами (numpy.random и scipy.stats)¶
Полное описание для scipy.stats доступно по ссылке. По ссылке можно прочитать полную документацию по работе с непрерывными (Continuous), дискретными (Discrete) и многомерными (Multivariate) распределениями. Пакет также предоставляет некоторое количество статистических методов, которые рассматриваются в курсах статистики. А описание для numpy.random находится по следующей ссылке.
import scipy.stats as sps
import numpy as np
import ipywidgets as widgets
import matplotlib.pyplot as plt
from numpy import random
import seaborn as sns
sns.set(palette='Set2', style='whitegrid', font_scale=1.3)
1. Работа с библиотекой numpy.random.¶
Базовые функции
randint(low[, high, size, dtype])— случайное число в диапазоне;choice(a[, size, replace, p])— Генерирует случайную выборку из заданного одномерного массива;shuffle(x)— Премешивает выходные значения массива;permutation(x)— Перемешивает значения самого массива и возвращает его перемешанным.
Создание простого случайного массива¶
Сгенерируем 5 случайных чисел из равномерного распределения [0, 1):
np.random.rand(5)
Создание случайных целых чисел:¶
Создадим случайный целочисленный массив размером 5x5 в диапазоне от 10 (включительно) до 20 (включительно)
np.random.randint(10, 20, (5, 5))
Также возможна генерация случайных чисел из конкретных распределений:¶
beta(a, b[, size])— бета-распределение;binomial(n, p[, size])— биномиальное распределение;exponential([scale, size])— экспоненциальное распределение;gamma(shape[, scale, size])— гамма-распределение;geometric(p[, size])— геометрическое распределение;normal([loc, scale, size])— нормальное распределение;poisson([lam, size])— пуассоновское распределение;uniform([low, high, size])— равномерное распределение.
Сгенерируем 5 случайных чисел из стандартного нормального распределения
np.random.normal(size=5)
Укажем параметры распределения
np.random.normal(70, 10, 5)
Генерация из пуассоновского распределения
np.random.poisson(lam=2.5, size=5)
Многомерное нормальное распределение:¶
Генерация случайных выборок из многомерного нормального распределения выполняется с помощью функции multivariate_normal(mean, cov[, size, ...]).
Сгенерируем 800 векторов из двумерного нормального распределения со средним [0, 0] и ковариационной матрицей [[6, -3], [-3, 3.5]].
mean = (1, 2)
cov = [[1, 0], [0, 1]]
x = np.random.multivariate_normal(mean, cov, (3, 3))
cov = np.array([[6, -3], [-3, 3.5]])
pts = np.random.multivariate_normal([0, 0], cov, size=800)
Мы можем визуализировать эти данные с помощью точечной диаграммы. Ориентация облака точек иллюстрирует отрицательную корреляцию компонентов этой выборки.
plt.figure(figsize=(6, 6))
plt.plot(pts[:, 0], pts[:, 1], '.', alpha = 0.5)
plt.axis('equal')
plt.show()
2. Работа с библиотекой scipy.stats.¶
Общий принцип:
Пусть $X$ — класс, реализующий некоторое распределение. Конкретное распределение с параметрами params можно получить как X(params). У него доступны следующие методы:
X(params).rvs(size=N)— генерация выборки размера $N$ (Random VariateS). Возвращаетnumpy.array;X(params).cdf(x)— значение функции распределения в точке $x$ (Cumulative Distribution Function);X(params).logcdf(x)— значение логарифма функции распределения в точке $x$;X(params).ppf(q)— $q$-квантиль (Percent Point Function);X(params).mean()— математическое ожидание;X(params).median()— медиана ($1/2$-квантиль);X(params).var()— дисперсия (Variance);X(params).std()— стандартное отклонение = корень из дисперсии (Standard Deviation).
Кроме того для непрерывных распределений определены функции
X(params).pdf(x)— значение плотности в точке $x$ (Probability Density Function);X(params).logpdf(x)— значение логарифма плотности в точке $x$.
А для дискретных
X(params).pmf(k)— значение дискретной плотности в точке $k$ (Probability Mass Function);X(params).logpdf(k)— значение логарифма дискретной плотности в точке $k$.
Все перечисленные выше методы применимы как к конкретному распределению X(params), так и к самому классу X. Во втором случае параметры передаются в сам метод. Например, вызов X.rvs(size=N, params) эквивалентен X(params).rvs(size=N).
Параметры могут быть следующими:
loc— параметр сдвига;scale— параметр масштаба;- и другие параметры (например, $n$ и $p$ для биномиального).
Для примера сгенерируем выборку размера $N = 200$ из распределения $\mathcal{N}(1, 9)$ и посчитаем некоторые статистики.
В терминах выше описанных функций у нас $X$ = sps.norm, а params = (loc=1, scale=3).
Примечание. Выборка — набор независимых одинаково распределенных случайных величин. Часто в разговорной речи выборку отождествляют с ее реализацией — значения случайных величин из выборки при "выпавшем" элементарном исходе.
sample = sps.norm(loc=1, scale=3).rvs(size=200)
print('Первые 10 значений выборки:\n', sample[:10])
print('Выборочное среденее: %.3f' % sample.mean())
print('Выборочная дисперсия: %.3f' % sample.var())
Вероятностные характеристики
print('Плотность:\t\t', sps.norm(loc=1, scale=3).pdf([-1, 0, 1, 2, 3]))
print('Функция распределения:\t', sps.norm(loc=1, scale=3).cdf([-1, 0, 1, 2, 3]))
$p$-квантиль распределения с функцией распределения $F$ — это число $min\{x: F(x) \geqslant p\}$.
print('Квантили:', sps.norm(loc=1, scale=3).ppf([0.05, 0.1, 0.5, 0.9, 0.95]))
Cгенерируем выборку размера $N = 200$ из распределения $Bin(10, 0.6)$ и посчитаем некоторые статистики.
В терминах выше описанных функций у нас $X$ = sps.binom, а params = (n=10, p=0.6).
sample = sps.binom(n=10, p=0.6).rvs(size=200)
print('Первые 10 значений выборки:\n', sample[:10])
print('Выборочное среденее: %.3f' % sample.mean())
print('Выборочная дисперсия: %.3f' % sample.var())
print('Дискретная плотность:\t', sps.binom(n=10, p=0.6).pmf([-1, 0, 5, 5.5, 10]))
print('Функция распределения:\t', sps.binom(n=10, p=0.6).cdf([-1, 0, 5, 5.5, 10]))
print('Квантили:', sps.binom(n=10, p=0.6).ppf([0.05, 0.1, 0.5, 0.9, 0.95]))
Отдельно есть класс для многомерного нормального распределения. Для примера сгенерируем выборку размера $N=200$ из распределения $\mathcal{N} \left( \begin{pmatrix} 1 \\ 1 \end{pmatrix}, \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} \right)$.
sample = sps.multivariate_normal(
mean=[1, 1], cov=[[2, 1], [1, 2]]
).rvs(size=200)
print('Первые 10 значений выборки:\n', sample[:10])
print('Выборочное среденее:', sample.mean(axis=0))
print('Выборочная матрица ковариаций:\n', np.cov(sample.T))
Некоторая хитрость :)
sample = sps.norm(loc=np.arange(10), scale=0.1).rvs(size=10)
print(sample)
Бывает так, что надо сгенерировать выборку из распределения, которого нет в `scipy.stats`.
Для этого надо создать класс, который будет наследоваться от класса rv_continuous для непрерывных случайных величин и от класса rv_discrete для дискретных случайных величин.
Пример из документации.
Для примера сгенерируем выборку из распределения с плотностью $f(x) = \frac{4}{15} x^3 I\{x \in [1, 2] = [a, b]\}$.
class cubic_gen(sps.rv_continuous):
def _pdf(self, x):
return 4 * x ** 3 / 15
cubic = cubic_gen(a=1, b=2, name='cubic')
sample = cubic.rvs(size=200)
print('Первые 10 значений выборки:\n', sample[:10])
print('Выборочное среденее: %.3f' % sample.mean())
print('Выборочная дисперсия: %.3f' % sample.var())
Если дискретная случайная величина может принимать небольшое число значений, то можно не создавать новый класс, как показано выше, а явно указать эти значения и из вероятности.
some_distribution = sps.rv_discrete(
name='some_distribution',
values=([1, 2, 3], [0.6, 0.1, 0.3]) # значения и вероятности
)
sample = some_distribution.rvs(size=200)
print('Первые 10 значений выборки:\n', sample[:10])
print('Выборочное среденее: %.3f' % sample.mean())
print('Частота значений по выборке:',
(sample == 1).mean(), (sample == 2).mean(), (sample == 3).mean())
3. Как генерируются случайные числа¶
Генерация случайных чисел разделена на два компонента: генератор битов и генератор случайных чисел.
BitGeneratorуправляет состоянием и предоставляет функции для создания случайных бинарных последовательностей.Генератор случайных чисел
Generatorпринимает набор из генератора битов, и преобразует их целевые распределения, например, нормальное распределение.
Генератор — пользовательский объект, который почти идентичен устаревшему RandomState . Он принимает экземпляр генератора битов в качестве аргумента. В настоящее время по умолчанию используется PCG64.
Можно воспользоваться функцие default_rng
rng = np.random.default_rng(12345)
print(rng)
print(rng.random())
Можно также создать экземпляр генератора напрямую с помощью экземпляра BitGenerator.
Чтобы использовать генератор бит PCG64 по умолчанию, можно создать его экземпляр напрямую и передать его генератору:
rng = np.random.Generator(np.random.PCG64(12345))
print(rng)
Аналогично использованию более старого генератора бит MT19937 (не рекомендуется), можно создать его экземпляр напрямую и передать его генератору:
rng = np.random.Generator(np.random.MT19937(12345))
print(rng)
Random.seed:¶
Данная функция используется для генерации одних и тех же случайных чисел при многократном выполнении кода на одной или разных машинах. По умолчанию генератор случайных чисел использует текущее системное время.
random.seed(svalue, version)
Параметр svalue является необязательным, и это начальное значение, необходимое для генерации случайного числа. Значение по умолчанию — None, и в таком случае генератор использует текущее системное время.
Если вы дважды используете одно и то же начальное значение, вы получите один и тот же результат
np.random.seed(10)
print(random.random())
np.random.seed(10)
print(random.random())
Другой пример, в котором мы генерируем одно и то же случайное число много раз
for i in range(5):
# устанавливаем опцию генератора
np.random.seed(11)
# генерируем число от 1 до 1000.
print(np.random.randint(1, 1000))
4. Производительность¶
Если вам нужно только сгенерировать случайный набор чисел без дополнительной работы с распределениями, рекомендуется использовать функции из numpy. Покажем это на примерах
size = 1000
Сравним время генерации 1000 случайных чисел двумя библиотеками
%timeit x = np.random.uniform(size=size, high=2)
%timeit x = sps.uniform.rvs(size=size, scale=2)
Как видим, numpy генерирует в 5 раз быстрее.
Если же мы будем иначе передавать параметры распределения в scipy, то будет еще дольше
%timeit x = sps.uniform(scale=2).rvs(size=size)
Это объясняется тем, что код sps.uniform(scale=2) создает объект распределения, на что тратится достаточно много времени
%timeit x = sps.uniform(scale=2)
Посмотрим также примеры с нормальным распределением
%timeit x = np.random.normal(size=size, loc=2, scale=5)
%timeit x = sps.norm.rvs(size=size, loc=2, scale=5)
%timeit x = sps.norm(loc=2, scale=5).rvs(size=size)
5. Свойства абсолютно непрерывных распределений¶
Прежде чем исследовать свойства распределений, напишем вспомогательную функцию для отрисовки плотности распределения.
def show_pdf(pdf, xmin, xmax, ymax, grid_size, distr_name, **kwargs):
"""
Рисует график плотности непрерывного распределения
pdf - плотность
xmin, xmax - границы графика по оси x
ymax - граница графика по оси y
grid_size - размер сетки, по которой рисуется график
distr_name - название распределения
kwargs - параметры плотности
"""
grid = np.linspace(xmin, xmax, grid_size)
plt.figure(figsize=(10, 5))
plt.plot(grid, pdf(grid, **kwargs), lw=5)
plt.grid(ls=':')
plt.xlabel('Значение', fontsize=18)
plt.ylabel('Плотность', fontsize=18)
plt.xlim((xmin, xmax))
plt.ylim((None, ymax))
title = 'Плотность {}'.format(distr_name)
plt.title(title.format(**kwargs), fontsize=20)
plt.show()
2.1 Нормальное распределение¶
$\mathcal{N}(a, \sigma^2)$ — нормальное распределение.
Параметры в scipy.stats:
loc= $a$,scale= $\sigma$.
Свойства распределения:
- математическое ожидание: $a$,
- дисперсия: $\sigma^2$.
Посмотрим, как выглядит плотность нормального стандартного распределения $\mathcal{N}(0, 1)$:
show_pdf(
pdf=sps.norm.pdf, xmin=-3, xmax=3, ymax=0.5, grid_size=100,
distr_name=r'$N({loc}, {scale})$', loc=0, scale=1
)
Сгенерируем значения из нормального стандартного распределения и сравним гистограмму с плотностью:
sample = sps.norm.rvs(size=1000) # выборка размера 1000
grid = np.linspace(-3, 3, 1000) # сетка для построения графика
plt.figure(figsize=(12, 5))
plt.hist(sample, bins=30, density=True,
alpha=0.6, label='Гистограмма выборки')
plt.plot(grid, sps.norm.pdf(grid), color='red',
lw=5, label='Плотность случайной величины')
plt.title(r'Случайная величина $\xi \sim \mathcal{N}$(0, 1)', fontsize=20)
plt.legend(fontsize=14, loc=1)
plt.show()
Исследуем, как меняется плотность распределения в зависимости от параметров.
# создать виджет, но не отображать его
ip = widgets.interactive(show_pdf,
pdf=widgets.fixed(sps.norm.pdf),
grid_size=widgets.IntSlider(min=25, max=300, step=25, value=100),
xmin=widgets.FloatSlider(min=-10, max=0, step=0.1, value=-5),
xmax=widgets.FloatSlider(min=0, max=10, step=0.1, value=5),
ymax=widgets.FloatSlider(min=0, max=2, step=0.1, value=1),
loc = widgets.FloatSlider(min=-10, max=10, step=0.1, value=0),
scale = widgets.FloatSlider(min=0.01, max=2, step=0.01, value=1),
distr_name = r'$N$({loc}, {scale})'
);
# отображаем слайдеры группами
display(widgets.HBox(ip.children[:2]))
display(widgets.HBox(ip.children[2:4]))
display(widgets.HBox(ip.children[5:7]))
# отображаем вывод функции
display(ip.children[-1])
ip.update() # чтобы функция запустилась до первого изменения слайдеров
Показательный пример с разными значениями параметров распределения:
grid = np.linspace(-7, 7, 1000) # сетка для построения графика
loc_values = [0, 3, 0] # набор значений параметра a
sigma_values = [1, 1, 2] # набор значений параметра sigma
plt.figure(figsize=(12, 6))
for i, (a, sigma) in enumerate(zip(loc_values, sigma_values)):
plt.plot(grid, sps.norm(a, sigma).pdf(grid), lw=5,
label='$\mathcal{N}' + '({}, {})$'.format(a, sigma))
plt.legend(fontsize=16)
plt.title('Плотности нормального распределения', fontsize=20)
plt.xlabel('Значение', fontsize=18)
plt.ylabel('Плотность', fontsize=18)
plt.show()
Значения параметров определяют положение и форму кривой на графике распределения, каждой комбинации параметров соответствует уникальное распределение.
Для нормального распределения:
- Параметр $loc = a$ отвечает за смещение кривой вдоль $\mathcal{Ox}$, тем самым определяя положение вертикальной оси симметрии плотности распределения. Вероятность того, что значение случайной величины $x$ попадет в отрезок $\mathcal{[m; n]}$, равна площади участка, зажатого кривой плотности, $\mathcal{Ox}$ и вертикальными прямыми ${x = m}$, ${x = n}$. В точке $a$ значение плотности распределения наибольшее, соответственно вероятность того, что значение случайной величины, имеющей нормальное распределение, попадет в окрестность точки $а$ — наибольшая.
- параметр ${scale = \sigma}$ отвечает за смещение экстремума вдоль $\mathcal{Oy}$ и "прижимание" кривой к вертикальной прямой ${x = a}$, тем самым увеличивая площадь под кривой плотности в окрестности точки $а$. Другими словами, этот параметр отвечает за дисперсию — меру разброса значений случайной величины. При уменьшении параметра $\sigma$ увеличивается вероятность того, что нормально распределенная случайная величина будет равна $a$. Это соответствует мере разброса значений случайной величины относительно её математического ожидания, то есть дисперсии $\sigma^2$.
Проверим несколько полезных свойств нормального распределения.
Пусть $\xi_1 \sim \mathcal{N}(a_1, \sigma_1^2)$ и $\xi_2 \sim \mathcal{N}(a_2, \sigma_2^2)$ — независимые случайные величины. Тогда $\xi_3 = \xi_1 + \xi_1 \sim \mathcal{N}(a_1 + a_2, \sigma_1^2 + \sigma_2^2)$
sample1 = sps.norm(loc=-1, scale=3).rvs(size=1000)
sample2 = sps.norm(loc=1, scale=4).rvs(size=1000)
sample3 = sample1 + sample2
grid = np.linspace(-15, 15, 1000)
plt.figure(figsize=(12, 5))
plt.hist(sample3, density=True, bins=30, alpha=0.6,
label=r'Гистограмма суммы $\xi_3 = \xi_1 + \xi_1$')
plt.plot(grid, sps.norm(-1 + 1, np.sqrt(3*3 + 4*4)).pdf(grid),
color='red', lw=5, label=r'Плотность $\mathcal{N}(0, 25)$')
plt.title(
r'Распределение $\xi_3=\xi_1+\xi_1\sim\mathcal{N}(-1 + 1, 3^2 + 4^2)$ ',
fontsize=20
)
plt.xlabel('Значение', fontsize=17)
plt.ylabel('Плотность', fontsize=17)
plt.legend(fontsize=16)
plt.show()
Пусть $\xi$ из $\mathcal{N}(a, \sigma^2)$. Тогда $\xi_{new} = c\xi\sim\mathcal{N}(c a, c^2 \sigma^2)$
sample = sps.norm(loc=5, scale=3).rvs(size=1000)
grid = np.linspace(-5, 30, 1000)
c = 2
new_sample = c*sample
plt.figure(figsize=(12,5))
plt.hist(new_sample, density=True, bins=30, alpha=0.6,
label=r'Гистограмма $\xi_{new} = c \xi$')
plt.plot(grid, sps.norm(c*5, c*3).pdf(grid), color='red',
lw=5, label=r'Плотность $\mathcal{N}(10, 36)$')
plt.title(
r'Распределение $\xi_{new}=c \xi\sim\mathcal{N}(2\cdot5, 4\cdot9)$',
fontsize=20
)
plt.xlabel('Значение', fontsize=17)
plt.ylabel('Плотность', fontsize=17)
plt.legend(fontsize=16)
plt.show()
2.2 Равномерное распределение¶
${U}(a, b)$ — равномерное распределение.
Параметры в scipy.stats:
loc= $a$,scale= $b-a$.
Свойства распределения:
- математическое ожидание: $\frac{a+b}{2}$,
- дисперсия: $\frac{(b-a)^2}{12}$.
show_pdf(
pdf=sps.uniform.pdf, xmin=-0.5, xmax=1.5, ymax=2, grid_size=10000,
distr_name=r'$U(0, 1)$', loc=0, scale=1
)
grid = np.linspace(-3, 3, 10001) # сетка для построения графика
sample = sps.uniform.rvs(size=1000) # выборка размера 1000
plt.figure(figsize=(12, 5))
plt.hist(sample, bins=30, density=True, alpha=0.6,
label='Гистограмма случайной величины')
plt.plot(grid, sps.uniform.pdf(grid), color='red', lw=5,
label='Плотность случайной величины')
plt.title(r'Случайная величина $\xi\sim U(0, 1)$', fontsize=20)
plt.xlim(-0.5, 1.5)
plt.legend(fontsize=14, loc=1)
plt.show()
# создать виджет, но не отображать его
ip = widgets.interactive(
show_pdf,
pdf=widgets.fixed(sps.uniform.pdf),
grid_size=widgets.IntSlider(min=25, max=300, step=25, value=500),
xmin=widgets.FloatSlider(min=-10, max=0, step=0.1, value=-5),
xmax=widgets.FloatSlider(min=0, max=10, step=0.1, value=5),
ymax=widgets.FloatSlider(min=0, max=2, step=0.1, value=1.4),
loc=widgets.FloatSlider(min=-4, max=0, step=0.1, value=0),
scale=widgets.FloatSlider(min=0.01, max=4, step=0.01, value=1),
distr_name=r'$U$({loc}, {loc} + {scale})'
);
# отображаем слайдеры группами
display(widgets.HBox(ip.children[:2]))
display(widgets.HBox(ip.children[2:4]))
display(widgets.HBox(ip.children[5:7]))
# отображаем вывод функции
display(ip.children[-1])
ip.update() # чтобы функция запустилась до первого изменения слайдеров
grid = np.linspace(-3, 7, 1000) # сетка для построения графика
loc_values = [0, 0, 4] # набор значений параметра a
scale_values = [1, 2, 1] # набор значений параметра scale
plt.figure(figsize=(12, 5))
for i, (loc, scale) in enumerate(zip(loc_values, scale_values)):
plt.plot(grid, sps.uniform(loc, scale).pdf(grid), lw=5, alpha=0.7,
label='$U' + '({}, {})$'.format(loc, scale))
plt.legend(fontsize=16)
plt.title('Плотности равномерного распределения', fontsize=20)
plt.xlabel('Значение', fontsize=18)
plt.ylabel('Плотность', fontsize=18)
plt.show()
Для равномерного распределения:
- Параметр ${loc = a}$ определяет начало отрезка, на котором случайная величина равномерно распределена.
- Параметр ${scale = b-a}$ определяет длину отрезка, на котором задана случайная величина. Значение плотности распределения на данном отрезке убывает с ростом данного параметра, то есть с ростом длины этого отрезка. Чем меньше длина отрезка, тем больше значение плотности вероятности на отрезке.
