Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- Лабораторная №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()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement