SHARE
TWEET

matstat lab1

allekco Nov 11th, 2019 89 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. Лабораторная №1. Проверка статистических гипотез
  2. Вариант № 1
  3.  
  4. Распределения:
  5.  
  6. 𝑋1  ~ N(1, 2) (объём выборки  𝑛1  = 100)
  7.  
  8. 𝑋2  ~ R(-1, 1) (объём выборки  𝑛2  = 100)
  9.  
  10. %matplotlib inline
  11. import numpy as np
  12. from scipy import stats
  13. import matplotlib.pyplot as plt
  14. Генерация выборок
  15. x = np.random.<distribution>(...params,size), где
  16.  
  17. distribution - распределение;
  18. ...params - параметры распределения;
  19. size - размер выборки
  20. Доступные распределения:
  21.  
  22. нормальное: normal(m, 𝜎2)
  23. равномерное: uniform(a, b)
  24. хи-квадрат: chisquare(k)
  25. # Размеры выборок
  26. n1 = 100
  27. n2 = 100
  28. # Функции для получения выборок
  29. def gen_x1():
  30.     return np.random.normal(1, 2, n1)
  31. def gen_x2():
  32.     return np.random.uniform(-1, 1, n2)
  33. # Конкретные выборки
  34. x1, x2 = gen_x1(), gen_x2()
  35. 1. Выборочные характеристики
  36. Необходимо:
  37.  
  38. Описать распределения 𝑋1 и 𝑋2, найти их МО и дисперсию, указать объём выборок
  39.  
  40. Рассчитать выборочные характеристики: среднее, 𝑠, 𝑠2
  41. Рассчитать выборочные характеристики для совокупной выборки 𝑥1 + 𝑥2
  42. def print_sample_chars(sample):
  43.     print("Среднее {}, s={}, s^2={}".format(
  44.         np.mean(sample),
  45.         np.std(sample),
  46.         np.var(sample)
  47.     ))
  48. print_sample_chars(x1)
  49. print_sample_chars(x2)
  50. pooled = np.concatenate([x1, x2])
  51. print_sample_chars(pooled)
  52. Среднее 0.8911644727059012, s=2.011210059192242, s^2=4.044965902196061
  53. Среднее -0.01004487405239921, s=0.5546407669037442, s^2=0.30762638031157347
  54. Среднее 0.440559799326751, s=1.5425111710859496, s^2=2.379340712924948
  55. Указания:
  56.  
  57. np.mean - среднее значение
  58. np.std - 𝑠 - оценка с.к.о.
  59. np.var - 𝑠2 - оценка дисперсии
  60. 2. Однопараметрические критерии
  61. Необходимо:
  62.  
  63. Для СВ 𝑋1 сформулировать гипотезы 𝐻0, проверяемые следующими тестами:
  64.  
  65. z-test
  66. t-test
  67. 𝜒2−𝑡𝑒𝑠𝑡 (𝑚 известно)
  68. 𝜒2−𝑡𝑒𝑠𝑡 (𝑚 неизвестно)
  69. Для каждой гипотезы рассчитать выборочное значение статистики критерия, p-value, выбрать уровень значимости 𝛼 и рассчитать ошибку статистического решения.
  70.  
  71. def criterion_z_test(sample):
  72.     m0 = 1
  73.     sigma = 2
  74.     mean = np.mean(sample)
  75.     n = len(sample)
  76.     return (mean - m0) / sigma * np.sqrt(n)
  77. def criterion_chi2_m(sample):
  78.     m0 = 1
  79.     sigma = 2
  80.     n = len(sample)
  81.     sum = 0
  82.     for i in range(n1):
  83.         sum = (sample[i] - m0)**2 + sum
  84.     S0_2 = sum/n
  85.     return n*S0_2/(sigma**2)
  86.    
  87. def criterion_chi2_no_m(sample):
  88.     m = np.mean(sample)
  89.     sigma = 2
  90.     n = len(sample)
  91.     sum = 0
  92.     for i in range(n):
  93.         sum = (sample[i] - m)**2 + sum
  94.     S_2 = sum/(n-1)
  95.     return (n-1)*S_2/(sigma**2)
  96. def criterion_t_test(sample):
  97.     n = len(sample)
  98.     s = np.std(sample) # s - оценка с.к.о.
  99.     mean = np.mean(sample) # выборочное среднее
  100.     m0 = 1 # основная гипотеза: МО генеральной совокупности для x2 составляет m0
  101.     return (mean - m0) / s * np.sqrt(n)
  102. z_test_dist = stats.norm
  103. t_test_dist = stats.t
  104. chi2_test_dist = stats.chi2
  105. chi2_test_f_dist = stats.f
  106. t-test
  107.  
  108. criterion_value = criterion_t_test(x1) # значение статистики критерия для гипотезы H0: m = m0, сигма неизвестна
  109. alpha = 0.05 # задаёмся уровнем значимости
  110. student_quantile = t_test_dist.ppf(1 - alpha / 2, n1 - 1) # рассчитываем квантиль распределения Стьюдента
  111. critical_value = student_quantile # критическое значение статистики критерия
  112. print("Значение критерия t-test: {}, критическое значение: {}".format(abs(criterion_value), critical_value))
  113. is_h0_true = abs(criterion_value) < critical_value
  114. print("Гипотеза H0 принимается" if is_h0_true else "Гипотеза H0 отклоняется")
  115. Значение критерия t-test: 0.5411445055013807, критическое значение: 1.9842169515086827
  116. Гипотеза H0 принимается
  117. def left_p_value(dist, criterion_value, sample_size):
  118.     return dist.cdf(criterion_value, sample_size - 1)
  119. def right_p_value(dist, criterion_value, sample_size):
  120.     return 1 - dist.cdf(criterion_value, sample_size - 1)
  121. def two_sided_p_value(dist, criterion_value, sample_size):
  122.     left_p = left_p_value(dist, criterion_value, sample_size)
  123.     return 2 * min(left_p, 1 - left_p)
  124. print("Двустороннее p-value: {}".format( two_sided_p_value(t_test_dist, criterion_value, n1) ))
  125. Двустороннее p-value: 0.5896236155621466
  126. z-test
  127.  
  128. criterion_value = criterion_z_test(x1)
  129. critical_value = z_test_dist.ppf(1 - alpha / 2, 0, 1)
  130. print("Значение критерия z-test: {}, критическое значение: {}".format(abs(criterion_value), critical_value))
  131. is_h0_true = abs(criterion_value) < critical_value
  132. print("Гипотеза H0 принимается" if is_h0_true else "Гипотеза H0 отклоняется")
  133. Значение критерия z-test: 0.5441776364704942, критическое значение: 1.959963984540054
  134. Гипотеза H0 принимается
  135. def left_p_value(dist, criterion_value, m, sigma):
  136.     return dist.cdf(criterion_value, m, sigma)
  137. def right_p_value(dist, criterion_value, m, sigma):
  138.     return 1 - dist.cdf(criterion_value, m, sigma)
  139. def two_sided_p_value(dist, criterion_value, m, sigma):
  140.     left_p = left_p_value(dist, criterion_value, m, sigma)
  141.     return 2 * min(left_p, 1 - left_p)
  142. print("Двустороннее p-value: {}".format( two_sided_p_value(z_test_dist, criterion_value, 0, 1) ))
  143. Двустороннее p-value: 0.5863192395510193
  144. 𝜒2−𝑡𝑒𝑠𝑡 (𝑚 известно)
  145.  
  146. criterion_value = criterion_chi2_m(x1)
  147. alpha = 0.05
  148. critical_value1 = chi2_test_dist.ppf(alpha/2, n1)
  149. critical_value2 = chi2_test_dist.ppf(1 - alpha/2, n1)
  150. critical_value = student_quantile
  151. print("Значение критерия chi2-test: {}, критическое значение: {}".format(abs(criterion_value), critical_value))
  152. is_h0_true = abs(criterion_value) < critical_value
  153. print("Гипотеза H0 принимается" if is_h0_true else "Гипотеза H0 отклоняется")
  154. Значение критерия chi2-test: 101.42027685493619, критическое значение: 1.9842169515086827
  155. Гипотеза H0 отклоняется
  156. def left_p_value(dist, criterion_value, sample_size):
  157.     return dist.cdf(criterion_value, sample_size )
  158. def right_p_value(dist, criterion_value, sample_size):
  159.     return 1 - dist.cdf(criterion_value, sample_size)
  160. def two_sided_p_value(dist, criterion_value, sample_size):
  161.     left_p = left_p_value(dist, criterion_value, sample_size)
  162.     return 2 * min(left_p, 1 - left_p)
  163. print("Двустороннее p-value: {}".format( two_sided_p_value(chi2_test_dist, criterion_value, n1) ))
  164. Двустороннее p-value: 0.8830809579840857
  165. 𝜒2−𝑡𝑒𝑠𝑡 (𝑚 неизвестно)
  166.  
  167. criterion_value = criterion_chi2_no_m(x1)
  168. alpha = 0.05
  169. critical_value1 = chi2_test_dist.ppf(alpha/2, n1)
  170. critical_value2 = chi2_test_dist.ppf(1 - alpha/2, n1)
  171. print("Значение критерия chi2 no m - test: {}, критическое значение: {}".format(abs(criterion_value), critical_value))
  172. is_h0_true = abs(criterion_value) < critical_value
  173. print("Гипотеза H0 принимается" if is_h0_true else "Гипотеза H0 отклоняется")
  174. Значение критерия chi2 no m - test: 101.1241475549015, критическое значение: 1.9842169515086827
  175. Гипотеза H0 отклоняется
  176. def left_p_value(dist, criterion_value, sample_size):
  177.     return dist.cdf(criterion_value, sample_size - 1)
  178. def right_p_value(dist, criterion_value, sample_size):
  179.     return 1 - dist.cdf(criterion_value, sample_size - 1)
  180. def two_sided_p_value(dist, criterion_value, sample_size):
  181.     left_p = left_p_value(dist, criterion_value, sample_size)
  182.     return 2 * min(left_p, 1 - left_p)
  183. print("Двустороннее p-value: {}".format( two_sided_p_value(chi2_test_dist, criterion_value, n1) ))
  184. Двустороннее p-value: 0.8436655076200767
  185. 3. Критерии для двух выборок
  186. Необходимо:
  187.  
  188. Выполнить задания пункта 2 для СВ 𝑋1 и 𝑋2 и следующих тестов:
  189.  
  190. 2-sample t-test
  191. 2-sample F-test (m известно)
  192. 2-sample F-test (m неизвестно)
  193. def criterion_t2_test(sample1, sample2):
  194.     s1 = np.std(sample1)
  195.     s2 = np.std(sample2)
  196.     n1 = len(sample1)
  197.     n2 = len(sample2)
  198.     S = (n1 - 1) * s1 * s1 + (n2 - 1) * s2 * s2
  199.     S /= n1 + n2 - 2
  200.    
  201.     m1 = np.mean(sample1)
  202.     m2 = np.mean(sample2)
  203.     z = (m1 - m2) / S
  204.     z /= np.sqrt(1.0 / n1 + 1.0 / n2)
  205.     return z
  206. def criterion_f_test_m(sample1, sample2):
  207.     S01_2 = 0
  208.     S02_2 = 0
  209.     m = 1
  210.     for i in range(len(sample1)):
  211.         S01_2 = S01_2 + (sample1[i]-m)**2
  212.     S01_2 = S01_2/len(sample1)
  213.     for i in range(len(sample2)):
  214.         S02_2 = S02_2 + (sample2[i]-m)**2
  215.     S02_2 = S02_2/len(sample2)
  216.     return S01_2/S02_2
  217. def criterion_f_test_no_m(sample1, sample2):
  218.     S1_2 = 0
  219.     S2_2 = 0
  220.     for i in range(len(sample1)):
  221.         S1_2 = S1_2 + (sample1[i] - np.mean(x1))**2
  222.     S1_2 = S1_2/(len(sample1) - 1)
  223.     for i in range(len(sample2)):
  224.         S2_2 = S2_2 + (sample2[i] - np.mean(x2))**2
  225.     S2_2 = S2_2/(len(sample2) - 1)
  226.     return S1_2/S2_2
  227. 2-sample t-test
  228.  
  229. criterion_value = criterion_t2_test(x1, x2)
  230. critical_value = t_test_dist.ppf(1 - alpha / 2, n1 + n2 - 2)
  231. print("Значение критерия: {}, критическое значение: {}".format(abs(criterion_value), critical_value))
  232. is_h0_true = abs(criterion_value) < critical_value
  233. print("Гипотеза H0 принимается" if is_h0_true else "Гипотеза H0 отклоняется")
  234. def left_p_value(dist, criterion_value, n1 , n2):
  235.     return dist.cdf(criterion_value, n1 + n2 - 2)
  236. def right_p_value(dist, criterion_value, n1 , n2):
  237.     return 1 - dist.cdf(criterion_value, n1 + n2 - 2)
  238. def two_sided_p_value(dist, criterion_value, n1 , n2):
  239.     left_p = left_p_value(dist, criterion_value, n1 , n2)
  240.     return 2 * min(left_p, 1 - left_p)
  241. print("Двустороннее p-value: {}".format( two_sided_p_value(t_test_dist, criterion_value_t2, n1, n2) ))
  242. Значение критерия: 2.9281457991023077, критическое значение: 1.9720174778338955
  243. Гипотеза H0 отклоняется
  244. Двустороннее p-value: 2.486254613298655e-07
  245. 2-sample F-test (m известно)
  246.  
  247. criterion_value = criterion_f_test_m(x1, x2)
  248. critical_value1 = chi2_test_f_dist.ppf(alpha/2, n1, n2)
  249. critical_value2 = chi2_test_f_dist.ppf(1 - alpha/2, n1, n2)
  250. print("Значение критерия: {}, критическое значение1: {}, критическое значение2: {}".format(criterion_value, critical_value1, critical_value2))
  251. is_h0_true1 = critical_value1 < criterion_value
  252. is_h0_true2 = criterion_value < critical_value2
  253. print("Гипотеза H0 принимается" if (is_h0_true1 and is_h0_true2) else "Гипотеза H0 отклоняется")
  254. def left_p_value(dist, criterion_value, n1, n2):
  255.     return dist.cdf(criterion_value, n1, n2)
  256. def right_p_value(dist, criterion_value, n1, n2):
  257.     return 1 - dist.cdf(criterion_value, n1, n2)
  258. def two_sided_p_value(dist, criterion_value, n1, n2):
  259.     left_p = left_p_value(dist, criterion_value, n1, n2)
  260.     return 2 * min(left_p, 1 - left_p)
  261. print("Двустороннее p-value: {}".format(two_sided_p_value(chi2_test_f_dist, criterion_value, n1, n2)))
  262. Значение критерия: 3.0552485688329765, критическое значение1: 0.674194729559777, критическое значение2: 1.4832509898927289
  263. Гипотеза H0 отклоняется
  264. Двустороннее p-value: 5.470469099932984e-08
  265. 2-sample F-test (m неизвестно)
  266.  
  267. criterion_value = criterion_f_test_no_m(x1, x2)
  268. critical_value1 = chi2_test_f_dist.ppf(alpha/2, n1, n2)
  269. critical_value2 = chi2_test_f_dist.ppf(1-alpha/2, n1, n2)
  270. print("Значение критерия: {}, критическое значение1: {}, критическое значение2: {}".format(criterion_value, critical_value1, critical_value2))
  271. is_h0_true1 = critical_value1 < criterion_value
  272. is_h0_true2 = criterion_value < critical_value2
  273. print("Гипотеза H0 принимается" if (is_h0_true1 and is_h0_true2) else "Гипотеза H0 отклоняется")
  274. def left_p_value(dist, criterion_value, n1, n2):
  275.     return dist.cdf(criterion_value, n1, n2)
  276. def right_p_value(dist, criterion_value, n1, n2):
  277.     return 1 - dist.cdf(criterion_value, n1, n2)
  278. def two_sided_p_value(dist, criterion_value, n1, n2):
  279.     left_p = left_p_value(dist, criterion_value, n1, n2)
  280.     return 2 * min(left_p, 1 - left_p)
  281. print("Двустороннее p-value: {}".format(two_sided_p_value(chi2_test_f_dist, criterion_value, n1 - 1, n2 - 1)))
  282. Значение критерия: 13.148956529993287, критическое значение1: 0.674194729559777, критическое значение2: 1.4832509898927289
  283. Гипотеза H0 отклоняется
  284. Двустороннее p-value: 2.220446049250313e-16
  285. 4. Выборочные характеристики для Z и P-value
  286. Необходимо: Вычислить для выборки мощностью 𝑁 выборочные значения для случайных величин p-value и значения статистики критерия. Гипотеза 𝐻0 задаётся вариантом лабораторной работы.
  287.  
  288. N = 500
  289. criteria = []
  290. p_values = []
  291. criterion_dist = t_test_dist
  292. for i in range(0, N):
  293.     sample = gen_x1()
  294.    
  295.     criterion_value = criterion_t_test(sample)
  296.     criteria.append(criterion_value)
  297.    
  298.     p_value_left = criterion_dist.cdf(criterion_value, n1 - 1)
  299.     p_value_right = 1 - p_value_left
  300.     p_value_two_sided = 2 * min(p_value_left, p_value_right)
  301.     p_values.append(p_value_two_sided)
  302. print("Критерий:")
  303. print("среднее = {}, s^2 = {}, s = {}".format(np.mean(criteria), np.var(criteria), np.std(criteria)))
  304. print("P-value:")
  305. print("среднее = {}, s^2 = {}, s = {}".format(np.mean(p_values), np.var(p_values), np.std(p_values)))
  306. Критерий:
  307. среднее = 0.03699143478992882, s^2 = 0.9962381650759466, s = 0.998117310277678
  308. P-value:
  309. среднее = 0.5016847546722268, s^2 = 0.08224295321751363, s = 0.28678032222855465
  310. Графики
  311. Графики теоретических распределений и гистограммы эспериментальных значений для статистики критерия и p-value:
  312.  
  313. fig, ax = plt.subplots()
  314. # histogram
  315. n, bins, patches = ax.hist(criteria, 20, density=True)
  316. # plot
  317. x = np.linspace(criterion_dist.ppf(0.05, n2 - 1), criterion_dist.ppf(0.95, n2 - 1), 100)
  318. ax.plot(x, criterion_dist.pdf(x, n2 - 1), 'r-', lw=5, alpha=0.6, label='Theoretical')
  319. ax.set_title("Статистика критерия")
  320. ax.set_ylabel("Плотность вероятности")
  321. fig.tight_layout()
  322. plt.show()
  323.  
  324. fig, ax = plt.subplots()
  325. # histogram
  326. n, bins, patches = ax.hist(p_values, 20, density=True)
  327. # plot
  328. x = np.linspace(stats.uniform.ppf(0.01), stats.uniform.ppf(0.99), 100)
  329. ax.plot(x, stats.uniform.pdf(x), 'r-', lw=5, alpha=0.6, label='Theoretical')
  330. ax.set_title("P-values")
  331. ax.set_ylabel("Плотность вероятности")
  332. fig.tight_layout()
  333. plt.show()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top