pyth_10_prob
Заметка 10. Численная теория вероятности.
курса Математический практикум по Питону.
Шокуров Антон В.
shokurov.anton.v@yandex.ru
http://машинноезрение.рф
Версия 0.16

Анотация

Вводятся базовые элементы питона (Python версии 3.xx) на базе ключевых библиотек относящихся к анализу данных: построение графиков/гистограмм и численой статистике. Последнее, в частности, используется для ввода ключевых понятий из анализа данных: выборка, плотность распределения и подгонка модели под данные.

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

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

Численная теория вероятностей

Генерирование случайных чисел.

Существуют специальные (библиотечные) функции для создания случайных чисел. Библиотека Numpy уже загружена и переименована как np. Поэтому возможен код:

In [1]:
import numpy as np
from matplotlib import pyplot as plt
import scipy.stats as models # Загружаем модуль знающий о сложных распределениях.
import scipy.special as sp
In [2]:
# В данном случае формируем числа согласно нормальному распределению.
np.random.randn( 4 ) # В скобках указывается размер массива.
Out[2]:
array([ 0.19568305,  1.24583303,  0.87614882, -0.49364795])
In [3]:
a = np.random.randn( 4, 3) # Можно размеры и по другим осям. В данном случае, будет матрица 4x3.
a
Out[3]:
array([[-0.13214215, -0.0944873 ,  0.66083757],
       [ 0.73859487,  0.58632269, -2.54962864],
       [ 1.0495145 ,  0.64112859, -1.2835252 ],
       [ 1.71876196, -1.54028453,  1.78735357]])
In [4]:
a.shape
Out[4]:
(4, 3)

np.random является модулем, т.е. он содержит набор функций со схожим смыслом. Равномерное распределение задается просто функцией rand из модуля np.random

In [5]:
np.random.rand(4) # Равномерное распределение на интервале [0, 1)
Out[5]:
array([0.44152324, 0.61974084, 0.81965149, 0.1667565 ])

Известно что по равномерной случайной величине можно сгенерировать произвольную слуайную величину использую обратную функцию.

In [6]:
bern = (np.random.rand(1000) < 0.7) + 0 # Распределение буернули с парметром 0.7. Зачем добавлен 0?)
bern[10:30]
Out[6]:
array([1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1])
In [7]:
np.sum(bern == 0)/len(bern), np.sum(bern == 1)/len(bern) # Численно получили почти то что требуется.
Out[7]:
(0.305, 0.695)
In [8]:
np.mean(bern == 0), np.mean(bern == 1) # Можно сразу посчитать среднее.
Out[8]:
(0.305, 0.695)

Дискретные случайные величины

In [9]:
data = np.random.randint(-10, 10, 10) # Равномерное целочисленное распределение между -10 и 10 невключительно.
data
Out[9]:
array([-10,  -7,  -8,  -3,  -6,   3,  -6,  -9,   6,  -8])

Давай посчитаем количество раз которое то или иной число встречается. Это можно сделать по разному.

In [10]:
equal = data == data[1] # Да, да, помимо стандартных арифметических операций над массивами
equal #  можно и логические операции делать.
Out[10]:
array([False,  True, False, False, False, False, False, False, False,
       False])
In [11]:
# А потом посчитать поличество "истин".
cnt = np.sum( equal ) # Напомню, при преобразовании в целое истина превращается в 1, а ложь в 0.
cnt
Out[11]:
1
In [12]:
cnt/len( equal ) # Долю можно найти поделив на общее количество чисел.
Out[12]:
0.1

Можно конечно посчитать и более сложные условия (диапазоны значений).

In [13]:
cnt_n = np.sum( data < 0 ) # Количество отрицательных
cnt_p = np.sum( data > 0 ) # Количество положительных
cnt_e = np.sum( data == 0 ) # Количество равных нулю
f'Количество положительных, отрицательных и равных нулю элементов соответственно равно {cnt_p}, {cnt_n} и {cnt_e}'
Out[13]:
'Количество положительных, отрицательных и равных нулю элементов соответственно равно 2, 8 и 0'

Для большей статестической правдоподобности.

In [14]:
data = np.random.randint(-10, 10, 100)

cnt_n = np.sum( data < 0 ) # Количество отрицательных
cnt_p = np.sum( data > 0 ) # Количество положительных
cnt_e = np.sum( data == 0 ) # Количество равных нулю

print( ">0 : ", cnt_p, "   <0 : ", cnt_n, " ==0 : ", cnt_e)
>0 :  37    <0 :  59  ==0 :  4
In [15]:
100/20 # Приблизительно равно количеству равных 0.
Out[15]:
5.0

Гистограмма

Давай посчитаем количество раз которое кстречается каждое из чисел. По научному это называется гистограммой.

In [16]:
# Ламбда функция для подсчета частоту встречаемости элемента в массиве.
cnt = lambda data, e : np.mean( data == e ) # Выборочная функция.
In [17]:
cnt( data, 0), cnt_e/len(data) # Вроде работает как надо.
Out[17]:
(0.04, 0.04)

Теперь нужно прогнать все возможные числа.

In [18]:
elems = range(-10, 10, 1) # все возможные исходы для равномерного распределения на рассматриваемом интервале.
list( elems ) # Выведем их для полноты понимания.
Out[18]:
[-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
In [19]:
hist = map( lambda e: cnt( data, e), elems )
hist = list( hist ) # Получили искомый список частот встречаемости.
In [20]:
n = 4 # Число для которого выводится частота встречаемости.
f'Число {n} встречается в выборке {hist[4]*100}%.'
Out[20]:
'Число 4 встречается в выборке 6.0%.'

Его можно по разному выводить. Можно просто списком. Но это не наглядно.

Можно графиком.

In [21]:
plt.plot( hist )
Out[21]:
[<matplotlib.lines.Line2D at 0x7f7b62dc0128>]

Но обычно все-таки использую спецаильный вид графиков -- гистограммы (прямоугольники).

In [22]:
plt.bar( elems, hist )
Out[22]:
<BarContainer object of 20 artists>
In [23]:
plt.bar( elems, hist )
plt.xticks( [-10, -6, -4, 0, 4, 6, 10] ) # Укажем системе где пометить колонки.
Out[23]:
([<matplotlib.axis.XTick at 0x7f7b62c58358>,
  <matplotlib.axis.XTick at 0x7f7b62ccd780>,
  <matplotlib.axis.XTick at 0x7f7b62ccdcc0>,
  <matplotlib.axis.XTick at 0x7f7b62c184a8>,
  <matplotlib.axis.XTick at 0x7f7b62c18978>,
  <matplotlib.axis.XTick at 0x7f7b62c18e48>,
  <matplotlib.axis.XTick at 0x7f7b62c22358>],
 <a list of 7 Text xticklabel objects>)
In [24]:
plt.bar( elems, hist )
plt.xticks( [-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 9] ); # ; Дабы не было печати.
In [25]:
# В питоне, а точнее в модуле numpy, есть готовая функция для этих целей: histogram.
hh = np.histogram( data ) # В аргументе подается сам массив данных.
hh
Out[25]:
(array([14,  6, 15, 18,  6,  4, 10,  8, 10,  9]),
 array([-10. ,  -8.1,  -6.2,  -4.3,  -2.4,  -0.5,   1.4,   3.3,   5.2,
          7.1,   9. ]))
In [26]:
# Изучим переменую hh, которая соответствует гистограмме.
type( hh ), len( hh )
# Она является парой (tuple)...содержит два элемента.
Out[26]:
(tuple, 2)
In [27]:
type( hh[0] ), type( hh[1] )
# состоящей из двух массивов (numpy.ndarray).
Out[27]:
(numpy.ndarray, numpy.ndarray)
In [28]:
hh # Смотрим содержимое переменной hist.
Out[28]:
(array([14,  6, 15, 18,  6,  4, 10,  8, 10,  9]),
 array([-10. ,  -8.1,  -6.2,  -4.3,  -2.4,  -0.5,   1.4,   3.3,   5.2,
          7.1,   9. ]))
In [29]:
hh[0] # Содержит встречаемость чисел в наборе данных .
Out[29]:
array([14,  6, 15, 18,  6,  4, 10,  8, 10,  9])
In [30]:
hh[1] # Содержит границы ячеек для которых посчитано количество.
Out[30]:
array([-10. ,  -8.1,  -6.2,  -4.3,  -2.4,  -0.5,   1.4,   3.3,   5.2,
         7.1,   9. ])
In [31]:
hh[0].shape, hh[1].shape # По умолчанию будет сделано 10 якеек. В данном случае это нам не годится.
Out[31]:
((10,), (11,))

Сформируем свои границы ячеек.

In [32]:
edges = list( elems )
edges + [10] # Дабы добавить правую границу (невключительно) последнего элемента.
Out[32]:
[-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
In [33]:
# Заново вызовем функцию передва границы ячеек во втором аргументе. Правая граница не включена в дапазон.
hh = np.histogram( data, edges + [10] ) # Поэтому нужно добавить 10 как самую праую границу.
hh # В случае целочисленных данных факически при подсчете используется только левая граница.
Out[33]:
(array([ 6,  8,  4,  2,  6,  9, 10,  8,  3,  3,  4,  0,  6,  4,  6,  2,  5,
         5,  2,  7]),
 array([-10,  -9,  -8,  -7,  -6,  -5,  -4,  -3,  -2,  -1,   0,   1,   2,
          3,   4,   5,   6,   7,   8,   9,  10]))
In [34]:
n = 4 # По аналогии.
f'Число {n} встречается в выборке {hh[0][4]} раз.' # При выводе вывелась дробная часть.
Out[34]:
'Число 4 встречается в выборке 6 раз.'

Возникает вопрос о том как убедится, что у нас результат совпадает.

In [35]:
hist = np.array( hist ) # Для удобства превратим список в массив.
In [36]:
hh[0].shape, hist.shape # Сначала убедимся, что совпадают размеры.
Out[36]:
((20,), (20,))
In [37]:
hh[0] # При малом объеме данных можно сравнить вывев массивы.
Out[37]:
array([ 6,  8,  4,  2,  6,  9, 10,  8,  3,  3,  4,  0,  6,  4,  6,  2,  5,
        5,  2,  7])
In [38]:
hist * len(data) # Умножаем на количество чисел ввиду того, что мы подсчитывали долю, а не количество.
Out[38]:
array([ 6.,  8.,  4.,  2.,  6.,  9., 10.,  8.,  3.,  3.,  4.,  0.,  6.,
        4.,  6.,  2.,  5.,  5.,  2.,  7.])

Можно сравнить визуально выведя обе гистограммы как график.

In [39]:
# Наивный способ отрисовки двух графиков столбцов.
plt.bar( elems, hist * len(data) ) 
plt.bar( elems, hh[0] ) # Первый график был закрашен вторым.
Out[39]:
<BarContainer object of 20 artists>
In [40]:
# Чтобы нарисовать два графика с колонками нужно изменить параметры отрисовки.
# Стандартно колонки имеют ширину 0.8. Аргумент width изменяет ширину колонки. Выставим её в половину.
plt.bar( elems, hist * len(data), width = 0.4, align = 'edge' ) # align задает способ выравнивания.
plt.bar( elems, hh[0], width = -0.4, align = 'edge' ) # Положительный по левую сторону
plt.xticks( [-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 9] ); # Отрицательный по правую сторону от метки.
plt.legend(["Наш подход", "Батарейка"])
Out[40]:
<matplotlib.legend.Legend at 0x7f7b62a2fd68>

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

Можно по-элементно сравнить массивы на совпадение.

In [41]:
eq = hh[0] == hist*100
eq # Опять же визуально видно, что все вроде истинно...
Out[41]:
array([ True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True, False])
In [42]:
eq.all( ) # Проверяет истинные ли все элементы вместе взятые.
Out[42]:
False
In [43]:
w = np.where( eq == False ) # Строим массив мест для которых ложь.
w
Out[43]:
(array([19]),)
In [44]:
hh[0][w], (hist*100)[w]
Out[44]:
(array([7]), array([7.]))
In [45]:
myhist = (hist*100).round() # Округляем числа до целого числа.
In [46]:
eq = hh[0] == myhist
eq
Out[46]:
array([ True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True])
In [47]:
eq.all() # Теперь должна дать истину.
Out[47]:
True
In [48]:
(hh[0] - hist*100).std() # Но можно и так, вычислив среднекадратичное отклонение разности.
Out[48]:
1.9357399876532503e-16
In [49]:
# Можно минуя функцию histogram сразу нарисовать график.
hhh = plt.hist( data, edges + [10] ) # При этом параметры отрисовки будут коректными автоматически.
In [50]:
# Вызовем функцию, которая вычисляет только гистограммы, она её не рисует.
hhh
Out[50]:
(array([ 6.,  8.,  4.,  2.,  6.,  9., 10.,  8.,  3.,  3.,  4.,  0.,  6.,
         4.,  6.,  2.,  5.,  5.,  2.,  7.]),
 array([-10,  -9,  -8,  -7,  -6,  -5,  -4,  -3,  -2,  -1,   0,   1,   2,
          3,   4,   5,   6,   7,   8,   9,  10]),
 <a list of 20 Patch objects>)

Видим, что результат совпадает. Скоре всего код функции plt.hist вызывает функцию np.histogram.

Непрерывные случайные величины

Рассмотрим непрерывную случайную величину.

In [51]:
data = np.random.rand( 100 ) # Массив из 100 чисел равномерно распределеных на интервале [0 1).
In [52]:
np.min( data ), np.max( data ) # Убедимся в этом.
Out[52]:
(0.003720135898781063, 0.999266461646274)

В данном случае мы не сможем посчитать количества раз которое та или иная величина встречается. Для адекватного подсчета нужно разбить интервал на равные подинтервалы.

In [53]:
# Допустим мы хотим разбить интервал [0, 1) на десять равных подинтервалов.
d = np.linspace(0, 1, 11) # Тогда нам нужно на этом интервале расположить 11 точек.
d
Out[53]:
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])
In [54]:
d[2], d[3] # Пример интервала.
Out[54]:
(0.2, 0.30000000000000004)
In [55]:
get = d[2] <= data # Это даст те числа которые больше равны d[2], т.е. 0.2. gt -- greater or equal than.
In [56]:
lt = data < d[3] # А это те которые строго меньше 0.3. l -- less than.
In [57]:
get & lt # Объеденить их можно так.
Out[57]:
array([ True, False, False, False, False, False, False, False, False,
        True, False,  True, False, False, False, False, False, False,
       False,  True, False, False, False, False, False, False, False,
       False, False, False, False, False,  True, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False,  True, False,  True, False,
       False, False, False,  True, False, False, False, False, False,
       False, False, False, False, False, False,  True, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
        True])
In [58]:
# Напишем функцию аналогичную дискретному случаю.
cnt = lambda data, e : np.mean( (e[0]<= data) & (data < e[1]) ) # Выборочная функция.
In [59]:
cnt( data, (d[2], d[3]) ), cnt( data, (d[5], d[6]) ) # Проведем подсчет для пары итервалов.
Out[59]:
(0.1, 0.09)
In [60]:
interval = list(map( lambda i: (d[i], d[i+1]), range(10) ) ) # Создаем список кортежей границ интервалов.
interval # Содержит кортежи интервалов.
Out[60]:
[(0.0, 0.1),
 (0.1, 0.2),
 (0.2, 0.30000000000000004),
 (0.30000000000000004, 0.4),
 (0.4, 0.5),
 (0.5, 0.6000000000000001),
 (0.6000000000000001, 0.7000000000000001),
 (0.7000000000000001, 0.8),
 (0.8, 0.9),
 (0.9, 1.0)]
In [61]:
hist2 = map( lambda e: cnt( data, e), interval)
hist2 = list(hist2)
hist2
Out[61]:
[0.08, 0.11, 0.1, 0.11, 0.1, 0.09, 0.09, 0.14, 0.07, 0.11]
In [62]:
plt.bar( d[:-1], hist2 ) # Какой то странный вид гистограммы для равномерного распределения....
Out[62]:
<BarContainer object of 10 artists>
In [63]:
plt.bar( d[:-1]+0.05, hist2, width=0.09 ) # Изменена ширена колонки дабы они не накладывались друг на друга.
plt.xticks( d[:-1]+0.05 ); # 0.05 Дабы центр колонки был там где надо.
In [64]:
plt.hist( data ) # Опять же можно воспользоватся готовой функцией.
Out[64]:
(array([ 8., 11., 10., 11., 11.,  8.,  9., 14.,  7., 11.]),
 array([0.00372014, 0.10327477, 0.2028294 , 0.30238403, 0.40193867,
        0.5014933 , 0.60104793, 0.70060256, 0.8001572 , 0.89971183,
        0.99926646]),
 <a list of 10 Patch objects>)

Плотность распределения

Возвращаемся к вопросу о формировании произвольного распределения по равномерному.

In [65]:
x = np.random.rand(1000)
nn = sp.erfinv( 2*x - 1  ) * np.sqrt(2)
In [66]:
plt.hist( nn )
Out[66]:
(array([ 11.,  28., 105., 142., 236., 211., 163.,  79.,  18.,   7.]),
 array([-2.8445062 , -2.25253848, -1.66057076, -1.06860304, -0.47663532,
         0.11533239,  0.70730011,  1.29926783,  1.89123555,  2.48320327,
         3.07517099]),
 <a list of 10 Patch objects>)
In [67]:
gg = models.gaussian_kde( nn ) # Можно приблизить плотность распределения за счет ядерного подхода.
x = np.linspace(-3, 3, 100)
plt.plot( x, gg(x) )
Out[67]:
[<matplotlib.lines.Line2D at 0x7f7b627ec278>]
In [68]:
1/np.sqrt(2*np.pi) # Максимальное значение плотности.
Out[68]:
0.3989422804014327
In [69]:
imax = np.argmax( gg(x) )
f'индекс {imax} в массиве x. В точке {x[imax]:.2} достигается максимальное значение {gg( x[imax] )[0]:.2}'
Out[69]:
'индекс 50 в массиве x. В точке 0.03 достигается максимальное значение 0.39'
In [70]:
gg = models.gaussian_kde( nn, 0.1 ) # Можно приблизить плотность распределения за счет ядерного подхода.
x = np.linspace(-3, 3, 100)
plt.plot( x, gg(x) )
Out[70]:
[<matplotlib.lines.Line2D at 0x7f7b57ae13c8>]
In [71]:
gg = models.gaussian_kde( nn, 0.01 ) # Можно приблизить плотность распределения за счет ядерного подхода.
x = np.linspace(-3, 3, 100)
plt.plot( x, gg(x) )
Out[71]:
[<matplotlib.lines.Line2D at 0x7f7b57ac63c8>]
In [ ]: