Лектор: Илья Щуров
Семинаристы: Евгения Ческидова, Евгений Ковалев
Ассистенты: Константин Ваниев, Софья Дымченко
На этом семинаре мы:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.stats as sts
%matplotlib inline
plt.style.use('ggplot')
# чтобы заработало, почему-то нужно писать в другой ячейке, отдельно от импортирования библиотек
plt.rcParams["figure.figsize"] = (10,7)
Зададим нормально распределенную случайную величину с параметрами $\mu = 1$ и $\sigma = 0.5$:
mu = 1
sigma = 0.5
# loc - параметр среднего, scale - параметр среднеквадратичного отклонения
norm_rv = sts.norm(loc=mu, scale=sigma)
Построим график функции распределения:
np.linspace(-1, 3, 100)
x = np.linspace(-1, 3, 100)
norm_cdf = norm_rv.cdf(x)
#plt.figure(figsize=(10,7))
plt.plot(x, norm_cdf)
plt.title('Normal RV CDF')
plt.xlabel('$x$', fontsize=16)
plt.ylabel('$F(x)$', fontsize=16)
plt.show()
Давайте добавим еще распределений!
# равномерное непрерывное
a = -1
b = 3
# loc - параметр левой границы, scale - параметр масштаба (то есть правая граница - loc+scale)
uniform_rv = sts.uniform(loc=a, scale=b-a)
# биномиальное
n = 5
p = 0.5
binom_rv = sts.binom(n=n, p=p)
# Пуассона
mu = 5
poisson_rv = sts.poisson(mu=mu)
x = np.linspace(-1, 10, 10000)
norm_cdf = norm_rv.cdf(x)
uniform_cdf = uniform_rv.cdf(x)
binom_cdf = binom_rv.cdf(x)
poisson_cdf = poisson_rv.cdf(x)
plt.plot(x, norm_cdf, label='normal')
plt.plot(x, uniform_cdf, label='uniform')
plt.plot(x, binom_cdf, label='binomial')
plt.plot(x, poisson_cdf, label='poisson')
plt.title('Multiple RV CDFs')
plt.xlabel('$x$', fontsize=16)
plt.ylabel('$F(x)$', fontsize=16)
plt.legend(loc='best')
plt.show()
x = np.linspace(-1, 10, 10000)
norm_pdf = norm_rv.pdf(x)
uniform_pdf = uniform_rv.pdf(x)
binom_pdf = binom_rv.pmf(x)
poisson_pdf = poisson_rv.pmf(x)
plt.plot(x, norm_pdf, label='normal')
plt.plot(x, uniform_pdf, label='uniform')
plt.plot(x, binom_pdf, label='binomial')
plt.plot(x, poisson_pdf, label='poisson')
plt.title('Multiple RV PDFs & PMFs')
plt.xlabel('$x$', fontsize=16)
plt.ylabel('$p(x)$', fontsize=16)
plt.legend(loc='best')
plt.show()
Сгенерируем выборку нормальной случайной величины, упомянутой ранее:
# выбираем число, задающее генератор случайных чисел, чтобы всегда получать один и тот же результат
np.random.seed(0)
norm_sample = norm_rv.rvs(size=1000)
Посчитайте с помощью этой выборки доверительный интервал для среднего с уровнем доверия $95\%$. Исходя из теоретических соображений, доверительный интервал с уровнем доверия $1 - \alpha$ может быть представлен в следующем виде: $$ \left(\hat{\mu} - z_{1 - \frac{\alpha}{2}}\frac{\hat{se}\left(\hat{\mu}\right)}{\sqrt{n}}, \ \hat{\mu} + z_{1 - \frac{\alpha}{2}}\frac{\hat{se}\left(\hat{\mu}\right)}{\sqrt{n}}\right) $$ Оценку $\hat{\mu}$ можно получить методом максимального правдоподобия - она будет равна выборочному среднему. Аналогично, оценка стандартной ошибки $\hat{se}\left(\hat{\mu}\right)$ будет выборочным стандартным отклонением. Квантиль же в нашем случае равен $z_{0.975}$.
Квантиль можно посчитать с помощью функции scipy.stats.norm.ppf()
def CI_calculate(sample):
# YOUR CODE HERE
return (CI_left_bound, CI_right_bound)
Содержит ли построенный интервал искомое значение среднего? А если провести эксперимент много раз, какая доля построенных доверительных интервалов будет содержать это значение? Убедитесь в этом на практике.
# YOUR CODE HERE
(btw, прикольная визуализация: http://rpsychologist.com/d3/CI/)
Проделайте похожий опыт, но с одновыборочным t-критерием Стьюдента. Сгенерируйте выборку той же нормальной случайной величины и проверьте гипотезу о равенстве средних.
Таблица критических значений t-критерия Стьюдента: http://www.statisticshowto.com/tables/t-distribution-table/
# YOUR CODE HERE
А что будет, если повторить эксперимент много раз? Какова доля (ошибочно) отвергнутых нулевых гипотез?
# YOUR CODE HERE
Проверим работу центральной предельной теоремы. Выберем какое-нибудь распределение, не очень похожее на нормальное - например, Распределение хи-квадрат. Построим случайную величину, имеющую данное распределение, с заданными параметрами:
# число степеней свободы
df = 4
chi2_rv = sts.chi2(df=df)
Сгенерируем выборку размера 1000:
sample = chi2_rv.rvs(size=1000)
sample[:10]
Построим гистограмму выборки, изобразив поверх нее теоретическую плотность распределения рассматриваемой случайной величины:
# сначала строим гистограмму выборки
plt.hist(sample, bins=20, density=True, label='sample points')
# теперь строим график плотности
x = np.linspace(0, 15, 10000)
chi2_pdf = chi2_rv.pdf(x)
plt.plot(x, chi2_pdf, label='density')
plt.xlabel('$x$')
plt.ylabel('number of sample points; $p(x)$')
plt.legend(loc='best')
plt.show()
Согласно ЦПТ, выборочное среднее выборки объема $n$ будет иметь распределение, близкое к $N\left(\mu, \frac{\sigma^2}{n}\right)$. Мы рассматриваем распределение хи-квадрат с $k=4$ степенями свободы. В этом случае $\mu = k = 4$, $\sigma^2 = 2k = 8$. Таким образом, $\bar{X}_n \to N\left(4, \frac{8}{n}\right)$ по распределению.
Теперь проверим это на практике. Будем генерировать набор выборок размера $n$ (например, 1000 выборок размера $n$), считать набор их выборочных средних, строить их гистограммы и соотносить их с теоретической плотностью соответствующего нормального распределения. Для этого напишем функцию:
def sample_mean_distr(n):
# YOUR CODE HERE
plt.show()
Посмотрим, что получается при подстановке различных значений $n$.
sample_mean_distr(5)
sample_mean_distr(10)
sample_mean_distr(25)
sample_mean_distr(50)
sample_mean_distr(100)
Какой вывод можно сделать из этого небольшого проделанного исследования?
На прошлом семинаре мы сделали предварительное предположение о том, что употребление алкоголя влияет на успехи в учебе: те, кто употребляет больше алкоголя в будние дни, учатся хуже, чем те, кто лучше себя в этом плане контролирует. Проверьте это предположение с помощью t-критерия Стьюдента.
data = pd.read_csv('math_students.csv', delimiter=',')
# YOUR CODE HERE
Рассмотрим четыре различных набора пар $\left(x_n, y_n\right)$.
data = pd.read_csv('anscombe.csv', index_col=0)
data
Посчитайте выборочное среднее и выборочную дисперсию для каждого столбца, не пользуясь функциями .mean(), .std() и .var() (и проверьте свой ответ с помощью них).
for col in data.columns:
# YOUR CODE HERE
print(col, col_mean, col_var)
Найдите коэффициент корреляции Пирсона для каждой пары $(x_n, y_n)$, не пользуясь функциями .corr(), .pearsonr() и тому подобными (и проверьте свой ответ с помощью них).
for i in range(4):
# YOUR CODE HERE
print(x_name, y_name, pearson_corr)
Теперь изобразите диаграмму рассеяния для каждой из пар $\left(x_n, y_n\right)$ вместе с прямой линейной регрессии $y = 3 + 0.5x$, обозначенной другим цветом.
for i in range(4):
# YOUR CODE HERE
plt.title('scatterplot #' + str(i + 1))
plt.xlabel(x_name)
plt.ylabel(y_name)
plt.show()