pyth_11_xi_statistics
Заметка 11. Статистики и Хи-квадрат.
курса Математический практикум по Питону.
Шокуров Антон В.
shokurov.anton.v@yandex.ru
http://машинноезрение.рф
Версия 0.17

Анотация

Исследуется использование компьютера для вычисления статистических параметров, например, вероятности события. Также рассматриваются статистики, которые используется при проверки гипотез.

Это предварительная версия! Любые замечания приветсвуются.

Слуйчайные величины и использование компьютера по назначению

Численная статистика

In [1]:
import numpy as np
from matplotlib import pyplot as plt
import scipy.stats as models # Загружаем модуль знающий о сложных распределениях.
import scipy.special as sp

Формирование случайной величины

In [2]:
gen_data = np.random.randn( 10000 ) * 3.3 + 0.5 # Вывода нет потому что мы присвоили результат переменой.
gen_data.shape
Out[2]:
(10000,)
In [3]:
# Строим гистограмму по заданному массиву чисел.
hist = plt.hist( gen_data, 50); # В переменную hist сохраняются сами значения гистограммы.
In [4]:
x = np.linspace(-10,10, 100)
p = models.norm.pdf( x, scale=3.3, loc=0.5 )
plt.plot(x, p)
Out[4]:
[<matplotlib.lines.Line2D at 0x7f753a31d0f0>]
In [5]:
plt.hist( gen_data, 50, density=True)
plt.plot( x, p, 'r')
Out[5]:
[<matplotlib.lines.Line2D at 0x7f753a2c24e0>]

Восстановление/оценка параметров случайной величины по выборке

In [6]:
mm, ss = models.norm.fit( gen_data ) # Вычислим параметры распределения по выборке.
mm, ss # Корректно когда искомое распределение известно.
Out[6]:
(0.5113507179058694, 3.285544200771083)
In [7]:
np.mean( gen_data ), np.std( gen_data )
Out[7]:
(0.5113507179058694, 3.285544200771083)
In [8]:
# Зная параметры распределения можно вычислить всякие вещи.
models.norm.cdf( -0.3, scale=ss, loc=mm ) # Например, можно вычислить вероятность события: < -0.3.
Out[8]:
0.40247515861613337

Мат ожидание от произвольной функции

In [9]:
def moment(x, n): # Вычислим моменты вручную.
    return x**n
In [10]:
models.norm.expect( lambda x:moment(x,1), scale=3.3, loc=0.5) # Через нашу функцию.
Out[10]:
0.49999999999999994
In [11]:
models.norm.moment( 1, scale=3.3, loc=0.5), models.norm.mean( scale=3.3, loc=0.5) # Аналитически.
Out[11]:
(0.5, 0.5)
In [12]:
models.norm.expect( lambda x:moment(x,2), scale=3.3, loc=0.5) # Через явную функцию.
Out[12]:
11.14
In [13]:
models.norm.moment( 2, scale=3.3, loc=0.5) # Аналитически, т.е. по формуле.
Out[13]:
11.139999999999999

Хм.. Какие же точки используются для интегрирования...

In [14]:
l = []
def dd(x):
    l.append(x)
    return x
In [15]:
models.norm.expect( dd, scale=3.3, loc=0.5)
Out[15]:
0.49999999999999994
In [16]:
plt.plot(l, np.zeros_like(l), '*')
Out[16]:
[<matplotlib.lines.Line2D at 0x7f753a1e4d68>]
In [17]:
models.norm.stats( scale=3.3, loc=0.5 ) # По формуле.
Out[17]:
(array(0.5), array(10.89))

Сравнение с численным

In [18]:
np.mean( gen_data < -0.3 ) # Если же нам неизвестно распределедние, то можно по данным.
Out[18]:
0.4029
In [19]:
models.norm.cdf( -0.3, scale=3.3, loc=0.5 ) # Если мы знаем случайную величину
Out[19]:
0.4042257258869503
In [20]:
gen_data = np.random.randn( 10000 ) * 3.3 + 0.5 # Вывода нет потому что мы присвоили результат переменой.
gen_data.shape
Out[20]:
(10000,)

a/b, t-statistics

Вычисления вдоль оси

In [21]:
a = np.array( [[2, 3], [-3, 4], [3, -2]]) # Введем матрицу.
a
Out[21]:
array([[ 2,  3],
       [-3,  4],
       [ 3, -2]])
In [22]:
# В Numpy есть ряд функций, которые выполняют операцию reduce.
np.sum( a ) # Например, вычисляет сумму массива a: 3 + -6 + -1 + 2.
Out[22]:
7

Такой вариант вызова выполняет операцию для всех элементов массива как для единого массива, т.е. без учета размерности. Но можно и указать вдоль какой именно размерности делать вычисления.

In [23]:
a[0,:], a[1,:], a[2,:]
Out[23]:
(array([2, 3]), array([-3,  4]), array([ 3, -2]))
In [24]:
a.sum( axis=1 )
Out[24]:
array([5, 1, 1])
In [25]:
a[:, 0], a[:, 1], a[:, 2] # Будет ошибка. Матрица то не квадратная.

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-25-ebf3a5647c89> in <module>
----> 1 a[:, 0], a[:, 1], a[:, 2] # Будет ошибка. Матрица то не квадратная.

IndexError: index 2 is out of bounds for axis 1 with size 2
In [26]:
a[:, 0], a[:, 1]
Out[26]:
(array([ 2, -3,  3]), array([ 3,  4, -2]))
In [27]:
a.sum( axis=0 )
Out[27]:
array([2, 5])

Для закрепления

In [28]:
a = np.array( [[[2, 3], [-3, 4]], [[-5, 7], [-9, 11]] ])
a.shape
Out[28]:
(2, 2, 2)
In [29]:
a[:, 1, 0]
Out[29]:
array([-3, -9])
In [30]:
np.sum( a, axis = 0) # Вдоль нулевой делаем и сохраняем результат в другие размерности ( матрица 2 на 2).
# Например, -3 + -9 равно 12 сохранено в элемент матрицы 1, 0.
Out[30]:
array([[ -3,  10],
       [-12,  15]])
In [31]:
a[0, :, 1]
Out[31]:
array([3, 4])
In [32]:
np.sum( a, axis = 1) # Суммируем вдоль первой и сохранаяем результат в оставшуюся размерность.
# Например, 3 + 4 равно 7 сохранено в элемент матрицы 0, 1.
Out[32]:
array([[ -1,   7],
       [-14,  18]])

Теперь продолжим наши изыскания...

In [33]:
n = 12
In [34]:
dat = np.random.randn(100000, n)*2 + 0.3
In [35]:
mm = dat.mean( axis=1 )
std = dat.std( axis=1 )
mm.shape, std.shape
Out[35]:
((100000,), (100000,))
In [36]:
def tstat_equal_var( means, stds, n ):
    s = np.sqrt( (stds[0]**2 + stds[1]**2)/n )
    v = (means[0] - means[1])/s
    return v
In [37]:
# Хотим получить доступ одновременно к двум выборкам (значениям их статистики).
mm = mm.reshape(-1,2) # Иначе прилось бы обращатся так: i*2:(i+1)*2
mm.shape
Out[37]:
(50000, 2)
In [38]:
mm[0]
Out[38]:
array([ 1.21026886, -0.41995588])
In [39]:
std = std.reshape(-1,2)
In [40]:
stats = list(zip( mm, std ))
stats[0]
Out[40]:
(array([ 1.21026886, -0.41995588]), array([2.65482984, 2.19860924]))
In [41]:
tt = list( map( lambda x: tstat_equal_var( x[0], x[1], n ) ,stats))
tt = np.array( tt )
In [42]:
plt.hist( tt, 30 );
In [43]:
np.mean( tt < -0.3 )
Out[43]:
0.39196
In [44]:
models.t.cdf( -0.3, 2 * n - 2 )
Out[44]:
0.38349673134391576
In [45]:
def gen_stat( dat ):
    mm = np.mean(dat, axis=1)
    std = np.std(dat, axis=1)
    mm = mm.reshape(-1,2)
    std = std.reshape(-1,2)
    stats = list(zip( mm, std ))
    tt = list( map( lambda x: tstat_equal_var( x[0], x[1], dat.shape[1] ) ,stats))
    return np.array( tt )
In [46]:
tt = gen_stat(np.random.randn(100000, n) + 5.3)
np.mean( tt < -0.3 )
Out[46]:
0.3915
In [47]:
tt = gen_stat(np.random.randn(100000, n)*3.5 - 5.5)
np.mean( tt < -0.3 )
Out[47]:
0.38432

Сравнение двух нормальных распределений с одинаковой дисперсией.

In [48]:
dat1 = np.random.randn(10000)*2 + 0.5
dat2 = np.random.randn(10000)*2 + 3.5
In [49]:
plt.hist( dat1, 50 );
plt.hist( dat2, 50 );
In [50]:
dat1 = np.random.randn(n)*2 + 0.5
dat2 = np.random.randn(n)*2 + 3.5
In [51]:
meann = np.mean( dat1 ), np.mean( dat2 )
stdd = np.std( dat1 ), np.std( dat2 )
In [52]:
v = tstat_equal_var( meann, stdd, n  )
v
Out[52]:
-7.0782527799604
In [53]:
vv = np.abs( v )
prob = np.mean( tt < -vv ) + np.mean( tt > vv )
prob
Out[53]:
0.0
In [54]:
models.t.ppf( prob/2,  2 * n - 2 )
Out[54]:
-inf
In [55]:
(1-models.t.cdf( vv, 2 * n - 2  ))
Out[55]:
2.1155557028418315e-07
In [56]:
models.t.ppf( 0.0032,  2 * n - 2 )
Out[56]:
-3.0129184824833284

Упражнение. (a/b test) Сделайте аналог для распределения Бернули. Рассматривайте, когда много зачений (считайте n большим). Тогда сработает закон больших чисел.

Хи квадрат -- Выражение из случайных величин

In [57]:
# Изучим Хи-квадрат численно.
gen_data = np.random.randn( 1000, 5) # 1000 раз по 5 нормальных распределений.
gen_data[0] # Так выглядит строчка матрицы.
Out[57]:
array([-0.58676315,  0.50620088,  0.03375566, -0.62145855,  0.74177084])
In [58]:
# Убедимся в правильности размера массива (матрицы).
gen_data.shape # В скобках перечислены размер массива вдоль каждого из измерений.
Out[58]:
(1000, 5)
In [59]:
gen_data[0].shape, gen_data[999].shape
Out[59]:
((5,), (5,))
In [60]:
sq_data = gen_data * gen_data # Массив состоящий их квадратов нормального распределения.
sq_data.shape
Out[60]:
(1000, 5)
In [61]:
gen_data2 = np.sum( sq_data, axis = 1 ) # сумируем вдоль оси с индекосм 1 (нумирация начинается с 0).
gen_data2.shape # получаем массив состоящий из суммы 5 квадратов номального распределения
# т.е. выборка из Хи-квадрат с параметром 5.
Out[61]:
(1000,)

Теперь забывае как мы его создали. Считаем, что данные откуда-то взялись.

In [62]:
hist = plt.hist( gen_data2, 50); # Сторим гистограмму. Она должна напоминать гистограмму Хи-квадрат.
In [63]:
cnt = np.mean( gen_data2 < 7  ) # Данная строчка посчитает сколько числе в выборке меньше 7.
cnt
Out[63]:
0.767
In [64]:
# Можно воспользоватся явной функцией (chi2), которая использует формулу.
models.chi2.cdf( 7, 5 ) # Первое число указывает точку (7), воторое значение параметра (5) распределения.
# Замечаем, что значение близко к рассчетному.
Out[64]:
0.7793596920632894
In [65]:
models.chi2.mean( 5 )
Out[65]:
5.0
In [66]:
it = models.chi2.interval( 0.95, 5 ) # Строим доверительный интервал.
it
Out[66]:
(0.831211613486663, 12.832501994030027)
In [67]:
models.chi2.cdf( it[0], 5 ), 1.0 - models.chi2.cdf( it[1], 5 ) # Проверяем доверительный интервал.
Out[67]:
(0.02500000000000003, 0.025000000000000022)
In [68]:
models.chi2.cdf( it[1], 5 ) - models.chi2.cdf( it[0], 5 ) # Проверяем его иначе.
Out[68]:
0.95
In [69]:
np.mean( (it[0] <= gen_data2) & (gen_data2 <= it[1]) ) # Проверяем численно.
Out[69]:
0.955

Критерий согласия Пирсона

In [70]:
N = 15
In [71]:
dat = models.poisson.rvs( mu=2.5, size=N )
dat
Out[71]:
array([3, 3, 2, 3, 2, 3, 2, 2, 3, 1, 1, 2, 1, 5, 3])
In [72]:
nmax = np.max( dat )
In [73]:
hh = np.histogram( dat, range(nmax+2), density=True )
hh
Out[73]:
(array([0.        , 0.2       , 0.33333333, 0.4       , 0.        ,
        0.06666667]), array([0, 1, 2, 3, 4, 5, 6]))
In [74]:
pp = models.poisson.pmf( mu=2.5, k=range(nmax+1) )
pp
Out[74]:
array([0.082085  , 0.2052125 , 0.25651562, 0.21376302, 0.13360189,
       0.06680094])
In [75]:
hh[0].shape, pp.shape
Out[75]:
((6,), (6,))
In [76]:
ss = np.sum( (hh[0] - pp)**2/pp )
ss
Out[76]:
0.40107929529730507
In [77]:
ss * N # Значение статистики Пирсона.
Out[77]:
6.016189429459576
In [78]:
# То с чем должны сравнить при 95% доверии. Если больше то отбрасываем гипотезу.
models.chi2.ppf( 0.95, N-1 )
Out[78]:
23.684791304840576
In [79]:
dat = models.poisson.rvs( mu=2.5, size=N )
nmax = np.max( dat )
hh = np.histogram( dat, range(nmax+2), density=True )
pp = models.poisson.pmf( mu=1.5, k=range(nmax+1) ) # Распределение с другим параметром.
ss = np.sum( (hh[0] - pp)**2/pp )
ss * N
Out[79]:
114.66520119386189
In [80]:
models.binom.rvs( p=0.3, n=3, size=N)
Out[80]:
array([1, 0, 1, 2, 1, 1, 2, 1, 1, 1, 1, 0, 1, 0, 0])
In [81]:
dat = models.randint.rvs( low=0, high=7, size=N) #binom.rvs( p=0.4, n=5, size=N)
nmax = np.max( dat )
hh = np.histogram( dat, range(nmax+2), density=True )
pp = models.poisson.pmf( mu=1.5, k=range(nmax+1) ) # Распределение с другим параметром.
ss = np.sum( (hh[0] - pp)**2/pp )
ss * N
Out[81]:
59.50018494425947
In [ ]: