matstat lab1

Лабораторная №1. Проверка статистических гипотез
Вариант № 1

Распределения:

𝑋1  ~ N(1, 2) (объём выборки  𝑛1  = 100)

𝑋2  ~ R(-1, 1) (объём выборки  𝑛2  = 100)

%matplotlib inline
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
Генерация выборок
x = np.random.<distribution>(...params,size), где

distribution - распределение;
...params - параметры распределения;
size - размер выборки
Доступные распределения:

нормальное: normal(m, 𝜎2)
равномерное: uniform(a, b)
хи-квадрат: chisquare(k)
# Размеры выборок
n1 = 100
n2 = 100
​
# Функции для получения выборок
def gen_x1():
    return np.random.normal(1, 2, n1)
​
def gen_x2():
    return np.random.uniform(-1, 1, n2)
​
# Конкретные выборки
x1, x2 = gen_x1(), gen_x2()
1. Выборочные характеристики
Необходимо:

Описать распределения 𝑋1 и 𝑋2, найти их МО и дисперсию, указать объём выборок

Рассчитать выборочные характеристики: среднее, 𝑠, 𝑠2
Рассчитать выборочные характеристики для совокупной выборки 𝑥1 + 𝑥2
def print_sample_chars(sample):
    print("Среднее {}, s={}, s^2={}".format(
        np.mean(sample),
        np.std(sample),
        np.var(sample)
    ))
​
print_sample_chars(x1)
print_sample_chars(x2)
pooled = np.concatenate([x1, x2])
print_sample_chars(pooled)
Среднее 0.8911644727059012, s=2.011210059192242, s^2=4.044965902196061
Среднее -0.01004487405239921, s=0.5546407669037442, s^2=0.30762638031157347
Среднее 0.440559799326751, s=1.5425111710859496, s^2=2.379340712924948
Указания:

np.mean - среднее значение
np.std - 𝑠 - оценка с.к.о.
np.var - 𝑠2 - оценка дисперсии
2. Однопараметрические критерии
Необходимо:

Для СВ 𝑋1 сформулировать гипотезы 𝐻0, проверяемые следующими тестами:

z-test
t-test
𝜒2−𝑡𝑒𝑠𝑡 (𝑚 известно)
𝜒2−𝑡𝑒𝑠𝑡 (𝑚 неизвестно)
Для каждой гипотезы рассчитать выборочное значение статистики критерия, p-value, выбрать уровень значимости 𝛼 и рассчитать ошибку статистического решения.

def criterion_z_test(sample):
    m0 = 1
    sigma = 2
    mean = np.mean(sample)
    n = len(sample)
    return (mean - m0) / sigma * np.sqrt(n)
​
def criterion_chi2_m(sample):
    m0 = 1
    sigma = 2
    n = len(sample)
    sum = 0
    for i in range(n1):
        sum = (sample[i] - m0)**2 + sum
    S0_2 = sum/n
    return n*S0_2/(sigma**2)

​
def criterion_chi2_no_m(sample):
    m = np.mean(sample)
    sigma = 2
    n = len(sample)
    sum = 0
    for i in range(n):
        sum = (sample[i] - m)**2 + sum
    S_2 = sum/(n-1)
    return (n-1)*S_2/(sigma**2)
​
def criterion_t_test(sample):
    n = len(sample)
    s = np.std(sample) # s - оценка с.к.о.
    mean = np.mean(sample) # выборочное среднее
    m0 = 1 # основная гипотеза: МО генеральной совокупности для x2 составляет m0
    return (mean - m0) / s * np.sqrt(n)
​
z_test_dist = stats.norm
t_test_dist = stats.t
chi2_test_dist = stats.chi2
chi2_test_f_dist = stats.f
​
t-test

criterion_value = criterion_t_test(x1) # значение статистики критерия для гипотезы H0: m = m0, сигма неизвестна
​
alpha = 0.05 # задаёмся уровнем значимости
student_quantile = t_test_dist.ppf(1 - alpha / 2, n1 - 1) # рассчитываем квантиль распределения Стьюдента
​
critical_value = student_quantile # критическое значение статистики критерия
​
print("Значение критерия t-test: {}, критическое значение: {}".format(abs(criterion_value), critical_value))
​
is_h0_true = abs(criterion_value) < critical_value
​
print("Гипотеза H0 принимается" if is_h0_true else "Гипотеза H0 отклоняется")
Значение критерия t-test: 0.5411445055013807, критическое значение: 1.9842169515086827
Гипотеза H0 принимается
def left_p_value(dist, criterion_value, sample_size):
    return dist.cdf(criterion_value, sample_size - 1)
​
def right_p_value(dist, criterion_value, sample_size):
    return 1 - dist.cdf(criterion_value, sample_size - 1)
​
def two_sided_p_value(dist, criterion_value, sample_size):
    left_p = left_p_value(dist, criterion_value, sample_size)
    return 2 * min(left_p, 1 - left_p)
​
print("Двустороннее p-value: {}".format( two_sided_p_value(t_test_dist, criterion_value, n1) ))
Двустороннее p-value: 0.5896236155621466
z-test

criterion_value = criterion_z_test(x1)
critical_value = z_test_dist.ppf(1 - alpha / 2, 0, 1)
​
print("Значение критерия z-test: {}, критическое значение: {}".format(abs(criterion_value), critical_value))
​
is_h0_true = abs(criterion_value) < critical_value
​
print("Гипотеза H0 принимается" if is_h0_true else "Гипотеза H0 отклоняется")
Значение критерия z-test: 0.5441776364704942, критическое значение: 1.959963984540054
Гипотеза H0 принимается
def left_p_value(dist, criterion_value, m, sigma):
    return dist.cdf(criterion_value, m, sigma)
​
def right_p_value(dist, criterion_value, m, sigma):
    return 1 - dist.cdf(criterion_value, m, sigma)
​
def two_sided_p_value(dist, criterion_value, m, sigma):
    left_p = left_p_value(dist, criterion_value, m, sigma)
    return 2 * min(left_p, 1 - left_p)
​
print("Двустороннее p-value: {}".format( two_sided_p_value(z_test_dist, criterion_value, 0, 1) ))
Двустороннее p-value: 0.5863192395510193
𝜒2−𝑡𝑒𝑠𝑡 (𝑚 известно)

criterion_value = criterion_chi2_m(x1)
alpha = 0.05
critical_value1 = chi2_test_dist.ppf(alpha/2, n1)
critical_value2 = chi2_test_dist.ppf(1 - alpha/2, n1)
critical_value = student_quantile
​
print("Значение критерия chi2-test: {}, критическое значение: {}".format(abs(criterion_value), critical_value))
​
is_h0_true = abs(criterion_value) < critical_value
​
print("Гипотеза H0 принимается" if is_h0_true else "Гипотеза H0 отклоняется")
Значение критерия chi2-test: 101.42027685493619, критическое значение: 1.9842169515086827
Гипотеза H0 отклоняется
def left_p_value(dist, criterion_value, sample_size):
    return dist.cdf(criterion_value, sample_size )
​
def right_p_value(dist, criterion_value, sample_size):
    return 1 - dist.cdf(criterion_value, sample_size)
​
def two_sided_p_value(dist, criterion_value, sample_size):
    left_p = left_p_value(dist, criterion_value, sample_size)
    return 2 * min(left_p, 1 - left_p)
​
print("Двустороннее p-value: {}".format( two_sided_p_value(chi2_test_dist, criterion_value, n1) ))
Двустороннее p-value: 0.8830809579840857
𝜒2−𝑡𝑒𝑠𝑡 (𝑚 неизвестно)

criterion_value = criterion_chi2_no_m(x1)
alpha = 0.05
critical_value1 = chi2_test_dist.ppf(alpha/2, n1)
critical_value2 = chi2_test_dist.ppf(1 - alpha/2, n1)
​
print("Значение критерия chi2 no m - test: {}, критическое значение: {}".format(abs(criterion_value), critical_value))
​
is_h0_true = abs(criterion_value) < critical_value
​
print("Гипотеза H0 принимается" if is_h0_true else "Гипотеза H0 отклоняется")
Значение критерия chi2 no m - test: 101.1241475549015, критическое значение: 1.9842169515086827
Гипотеза H0 отклоняется
def left_p_value(dist, criterion_value, sample_size):
    return dist.cdf(criterion_value, sample_size - 1)
​
def right_p_value(dist, criterion_value, sample_size):
    return 1 - dist.cdf(criterion_value, sample_size - 1)
​
def two_sided_p_value(dist, criterion_value, sample_size):
    left_p = left_p_value(dist, criterion_value, sample_size)
    return 2 * min(left_p, 1 - left_p)
​
print("Двустороннее p-value: {}".format( two_sided_p_value(chi2_test_dist, criterion_value, n1) ))
Двустороннее p-value: 0.8436655076200767
3. Критерии для двух выборок
Необходимо:

Выполнить задания пункта 2 для СВ 𝑋1 и 𝑋2 и следующих тестов:

2-sample t-test
2-sample F-test (m известно)
2-sample F-test (m неизвестно)
def criterion_t2_test(sample1, sample2):
    s1 = np.std(sample1)
    s2 = np.std(sample2)
    n1 = len(sample1)
    n2 = len(sample2)
    S = (n1 - 1) * s1 * s1 + (n2 - 1) * s2 * s2
    S /= n1 + n2 - 2

    m1 = np.mean(sample1)
    m2 = np.mean(sample2)
    z = (m1 - m2) / S
    z /= np.sqrt(1.0 / n1 + 1.0 / n2)
    return z
​
def criterion_f_test_m(sample1, sample2):
    S01_2 = 0
    S02_2 = 0
    m = 1
    for i in range(len(sample1)):
        S01_2 = S01_2 + (sample1[i]-m)**2
    S01_2 = S01_2/len(sample1)
    for i in range(len(sample2)):
        S02_2 = S02_2 + (sample2[i]-m)**2
    S02_2 = S02_2/len(sample2)
    return S01_2/S02_2
​
def criterion_f_test_no_m(sample1, sample2):
    S1_2 = 0
    S2_2 = 0
    for i in range(len(sample1)):
        S1_2 = S1_2 + (sample1[i] - np.mean(x1))**2
    S1_2 = S1_2/(len(sample1) - 1)
    for i in range(len(sample2)):
        S2_2 = S2_2 + (sample2[i] - np.mean(x2))**2
    S2_2 = S2_2/(len(sample2) - 1)
    return S1_2/S2_2
2-sample t-test

criterion_value = criterion_t2_test(x1, x2)
​
critical_value = t_test_dist.ppf(1 - alpha / 2, n1 + n2 - 2)
​
print("Значение критерия: {}, критическое значение: {}".format(abs(criterion_value), critical_value))
​
is_h0_true = abs(criterion_value) < critical_value
​
print("Гипотеза H0 принимается" if is_h0_true else "Гипотеза H0 отклоняется")
​
​
def left_p_value(dist, criterion_value, n1 , n2):
    return dist.cdf(criterion_value, n1 + n2 - 2)
​
def right_p_value(dist, criterion_value, n1 , n2):
    return 1 - dist.cdf(criterion_value, n1 + n2 - 2)
​
def two_sided_p_value(dist, criterion_value, n1 , n2):
    left_p = left_p_value(dist, criterion_value, n1 , n2)
    return 2 * min(left_p, 1 - left_p)
​
print("Двустороннее p-value: {}".format( two_sided_p_value(t_test_dist, criterion_value_t2, n1, n2) ))
Значение критерия: 2.9281457991023077, критическое значение: 1.9720174778338955
Гипотеза H0 отклоняется
Двустороннее p-value: 2.486254613298655e-07
2-sample F-test (m известно)

criterion_value = criterion_f_test_m(x1, x2)
​
critical_value1 = chi2_test_f_dist.ppf(alpha/2, n1, n2)
​
critical_value2 = chi2_test_f_dist.ppf(1 - alpha/2, n1, n2)
​
print("Значение критерия: {}, критическое значение1: {}, критическое значение2: {}".format(criterion_value, critical_value1, critical_value2))
​
is_h0_true1 = critical_value1 < criterion_value
is_h0_true2 = criterion_value < critical_value2
​
print("Гипотеза H0 принимается" if (is_h0_true1 and is_h0_true2) else "Гипотеза H0 отклоняется")
​
def left_p_value(dist, criterion_value, n1, n2):
    return dist.cdf(criterion_value, n1, n2)
​
def right_p_value(dist, criterion_value, n1, n2):
    return 1 - dist.cdf(criterion_value, n1, n2)
​
def two_sided_p_value(dist, criterion_value, n1, n2):
    left_p = left_p_value(dist, criterion_value, n1, n2)
    return 2 * min(left_p, 1 - left_p)
​
print("Двустороннее p-value: {}".format(two_sided_p_value(chi2_test_f_dist, criterion_value, n1, n2)))
​
Значение критерия: 3.0552485688329765, критическое значение1: 0.674194729559777, критическое значение2: 1.4832509898927289
Гипотеза H0 отклоняется
Двустороннее p-value: 5.470469099932984e-08
2-sample F-test (m неизвестно)

criterion_value = criterion_f_test_no_m(x1, x2)
​
critical_value1 = chi2_test_f_dist.ppf(alpha/2, n1, n2)
​
critical_value2 = chi2_test_f_dist.ppf(1-alpha/2, n1, n2)
​
print("Значение критерия: {}, критическое значение1: {}, критическое значение2: {}".format(criterion_value, critical_value1, critical_value2))
​
is_h0_true1 = critical_value1 < criterion_value
is_h0_true2 = criterion_value < critical_value2
​
print("Гипотеза H0 принимается" if (is_h0_true1 and is_h0_true2) else "Гипотеза H0 отклоняется")
​
def left_p_value(dist, criterion_value, n1, n2):
    return dist.cdf(criterion_value, n1, n2)
​
def right_p_value(dist, criterion_value, n1, n2):
    return 1 - dist.cdf(criterion_value, n1, n2)
​
def two_sided_p_value(dist, criterion_value, n1, n2):
    left_p = left_p_value(dist, criterion_value, n1, n2)
    return 2 * min(left_p, 1 - left_p)
​
print("Двустороннее p-value: {}".format(two_sided_p_value(chi2_test_f_dist, criterion_value, n1 - 1, n2 - 1)))
Значение критерия: 13.148956529993287, критическое значение1: 0.674194729559777, критическое значение2: 1.4832509898927289
Гипотеза H0 отклоняется
Двустороннее p-value: 2.220446049250313e-16
4. Выборочные характеристики для Z и P-value
Необходимо: Вычислить для выборки мощностью 𝑁 выборочные значения для случайных величин p-value и значения статистики критерия. Гипотеза 𝐻0 задаётся вариантом лабораторной работы.

N = 500
​
criteria = []
p_values = []
​
criterion_dist = t_test_dist
​
for i in range(0, N):
    sample = gen_x1()

    criterion_value = criterion_t_test(sample)
    criteria.append(criterion_value)

    p_value_left = criterion_dist.cdf(criterion_value, n1 - 1)
    p_value_right = 1 - p_value_left
    p_value_two_sided = 2 * min(p_value_left, p_value_right)
    p_values.append(p_value_two_sided)
​
print("Критерий:")
print("среднее = {}, s^2 = {}, s = {}".format(np.mean(criteria), np.var(criteria), np.std(criteria)))
​
print("P-value:")
print("среднее = {}, s^2 = {}, s = {}".format(np.mean(p_values), np.var(p_values), np.std(p_values)))
Критерий:
среднее = 0.03699143478992882, s^2 = 0.9962381650759466, s = 0.998117310277678
P-value:
среднее = 0.5016847546722268, s^2 = 0.08224295321751363, s = 0.28678032222855465
Графики
Графики теоретических распределений и гистограммы эспериментальных значений для статистики критерия и p-value:

fig, ax = plt.subplots()
​
# histogram
n, bins, patches = ax.hist(criteria, 20, density=True)
​
# plot
x = np.linspace(criterion_dist.ppf(0.05, n2 - 1), criterion_dist.ppf(0.95, n2 - 1), 100)
ax.plot(x, criterion_dist.pdf(x, n2 - 1), 'r-', lw=5, alpha=0.6, label='Theoretical')
​
ax.set_title("Статистика критерия")
ax.set_ylabel("Плотность вероятности")
​
fig.tight_layout()
plt.show()

fig, ax = plt.subplots()
​
# histogram
n, bins, patches = ax.hist(p_values, 20, density=True)
​
# plot
x = np.linspace(stats.uniform.ppf(0.01), stats.uniform.ppf(0.99), 100)
ax.plot(x, stats.uniform.pdf(x), 'r-', lw=5, alpha=0.6, label='Theoretical')
​
ax.set_title("P-values")
ax.set_ylabel("Плотность вероятности")
​
fig.tight_layout()
plt.show()