In [23]:
plt.plot( hist )
Out[23]:
[<matplotlib.lines.Line2D at 0x7fe9e499ee48>]
Математический практикум по Питону.
Вводятся базовые элементы питона (Python версии 3.xx) на базе ключевых библиотек относящихся к анализу данных: построение графиков/гистограмм и численной статистике. Последнее, в частности, используется для ввода ключевых понятий из анализа данных: выборка, плотность распределения и подгонка модели под данные.
Это предварительная версия! Любые замечания приветствуются.
import numpy as np
from matplotlib import pyplot as plt
Нормальное
Существуют специальные (библиотечные) функции для создания случайных чисел. Библиотека Numpy уже загружена и переименована как np. Поэтому возможен код:
# В данном случае формируем числа согласно нормальному распределению.
np.random.randn( 4 ) # В скобках указывается размер массива.
array([-0.36444611, 0.85788507, -1.39084404, -0.62661321])
a = np.random.randn( 4, 3) # Можно размеры и по другим осям. В данном случае, будет матрица 4x3.
a
array([[ 0.79334778, -0.86260596, 1.91965298],
[ 0.13963103, 2.55403883, 1.16350701],
[ 0.44533086, 1.15812127, -0.96379915],
[ 1.76622369, 0.56193136, -0.3312121 ]])
a.shape
(4, 3)
np.random является модулем, т.е. он содержит набор функций со схожим смыслом. Равномерное распределение задается просто функцией rand из модуля np.random
Равномерное
np.random.rand(4) # Равномерное распределение на интервале [0, 1)
array([0.9701337 , 0.83863067, 0.83016209, 0.71468217])
Известно что по равномерной случайной величине можно сгенерировать произвольную слуайную величину использую обратную функцию.
bern = (np.random.rand(1000) < 0.7) + 0 # Распределение буернули с парметром 0.7. Зачем добавлен 0?)
bern[10:30]
array([1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1])
np.sum(bern == 0)/len(bern), np.sum(bern == 1)/len(bern) # Численно получили почти то что требуется.
(0.294, 0.706)
np.mean(bern == 0), np.mean(bern == 1) # Можно сразу посчитать среднее.
(0.294, 0.706)
data = np.random.randint(-10, 10, 10) # Равномерное целочисленное распределение между -10 и 10 невключительно.
data
array([ 8, -5, 7, -10, 0, -6, -9, -3, -8, 3])
Давай посчитаем количество раз которое то или иной число встречается. Это можно сделать по разному.
equal = data == data[2] # Да, да, помимо стандартных арифметических операций над массивами
equal # можно и логические операции делать.
array([False, False, True, False, False, False, False, False, False,
False])
# А потом посчитать поличество "истин".
cnt = np.sum( equal ) # Напомню, при преобразовании в целое истина превращается в 1, а ложь в 0.
cnt
1
cnt/len( equal ) # Долю можно найти поделив на общее количество чисел.
0.1
Можно конечно посчитать и более сложные условия (диапазоны значений).
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}'
'Количество положительных, отрицательных и равных нулю элементов соответственно равно 3, 6 и 1'
Для большей статестической правдоподобности.
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 : 34 <0 : 57 ==0 : 9
100/20 # Приблизительно равно количеству равных 0.
5.0
np.max(data) - np.min(data)+1
20
Давай посчитаем количество раз которое кстречается каждое из чисел. По научному это называется гистограммой.
# Ламбда функция для подсчета частоту встречаемости элемента в массиве.
cnt = lambda data, e : np.mean( data == e ) # Выборочная функция.
cnt( data, 0), cnt_e/len(data) # Вроде работает как надо.
(0.09, 0.09)
cnt(data,5)
0.05
Теперь нужно прогнать все возможные числа.
# все возможные исходы для равномерного распределения
elems = range(-10, 10, 1) # на рассматриваемом интервале.
list( elems ) # Выведем их для полноты понимания.
[-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
hist = map( lambda e: cnt( data, e), elems )
hist = list( hist ) # Получили искомый список частот встречаемости.
n = 4 # Число для которого выводится частота встречаемости.
f'Число {n} встречается в выборке {hist[10+n]*100}%.'
'Число 4 встречается в выборке 4.0%.'
Его можно по разному выводить. Можно просто списком. Но это не наглядно.
Можно графиком.
plt.plot( hist )
[<matplotlib.lines.Line2D at 0x7fe9e499ee48>]
Но обычно все-таки использую спецаильный вид графиков -- гистограммы (прямоугольники).
plt.bar( elems, hist )
<BarContainer object of 20 artists>
z=np.linspace(-5,5, 20)
plt.bar(z, z*z)
<BarContainer object of 20 artists>
plt.bar( elems, hist )
plt.xticks( [-10, -9, -8, -6, -4, 0, 4, 6, 10] ); # Укажем системе где пометить колонки.
plt.bar( elems, hist )
plt.xticks( [-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 9] ); # ; Дабы не было печати.
plt.bar( elems, hist )
plt.xticks( elems ); # Укажем системе где пометить колонки.
# В питоне, а точнее в модуле numpy, есть готовая функция для этих целей: histogram.
hh = np.histogram( data ) # В аргументе подается сам массив данных.
hh
(array([14, 10, 10, 7, 16, 10, 12, 9, 6, 6]),
array([-10. , -8.1, -6.2, -4.3, -2.4, -0.5, 1.4, 3.3, 5.2,
7.1, 9. ]))
np.unique( data )#.shape
array([-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2,
3, 4, 5, 6, 7, 8, 9])
# Изучим переменую hh, которая соответствует гистограмме.
type( hh ), len( hh )
# Она является парой (tuple)...содержит два элемента.
(tuple, 2)
type( hh[0] ), type( hh[1] )
# состоящей из двух массивов (numpy.ndarray).
(numpy.ndarray, numpy.ndarray)
hh[0].shape, hh[1].shape
((10,), (11,))
hh[0].dtype, hh[1].dtype
(dtype('int64'), dtype('float64'))
hh # Смотрим содержимое переменной hist.
(array([14, 10, 10, 7, 16, 10, 12, 9, 6, 6]),
array([-10. , -8.1, -6.2, -4.3, -2.4, -0.5, 1.4, 3.3, 5.2,
7.1, 9. ]))
hh[0] # Содержит встречаемость чисел в наборе данных .
array([14, 10, 10, 7, 16, 10, 12, 9, 6, 6])
# Содержит границы ячеек для которых посчитано количество.
hh[1]
array([-10. , -8.1, -6.2, -4.3, -2.4, -0.5, 1.4, 3.3, 5.2,
7.1, 9. ])
hh[0].shape, hh[1].shape # По умолчанию будет сделано 10 якеек. В данном случае это нам не годится.
((10,), (11,))
Сформируем свои границы ячеек.
edges = list( elems )
edges + [10] # Дабы добавить правую границу (невключительно) последнего элемента.
[-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
# Заново вызовем функцию передва границы ячеек во втором аргументе. Правая граница не включена в дапазон.
hh = np.histogram( data, edges + [10] ) # Поэтому нужно добавить 10 как самую праую границу.
hh # В случае целочисленных данных факически при подсчете используется только левая граница.
(array([ 7, 7, 5, 5, 4, 6, 2, 5, 6, 10, 9, 1, 6, 6, 4, 5, 4,
2, 4, 2]),
array([-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2,
3, 4, 5, 6, 7, 8, 9, 10]))
n = 4 # По аналогии.
f'Число {n} встречается в выборке {hh[0][10+n]} раз.' # При выводе вывелась дробная часть.
'Число 4 встречается в выборке 4 раз.'
Построение гистограммы
# Можно минуя функцию histogram сразу нарисовать график.
hhh = plt.hist( data, edges + [10] ) # При этом параметры отрисовки будут коректными автоматически.
plt.xticks( edges + [10] );
# Вызовем функцию, которая вычисляет только гистограммы, она её не рисует.
hhh
(array([ 7., 7., 5., 5., 4., 6., 2., 5., 6., 10., 9., 1., 6.,
6., 4., 5., 4., 2., 4., 2.]),
array([-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2,
3, 4, 5, 6, 7, 8, 9, 10]),
<BarContainer object of 20 artists>)
Видим, что результат совпадает. Скоре всего код функции plt.hist вызывает функцию np.histogram.
Упр. Вычислить гистограмму совсем вручную: использую цикл и список...
Возникает вопрос о том как убедится, что у нас результат совпадает.
hist = np.array( hist ) # Для удобства превратим список в массив.
hh[0].shape, hist.shape # Сначала убедимся, что совпадают размеры.
((20,), (20,))
hh[0] # При малом объеме данных можно сравнить выведевв массивы.
array([ 7, 7, 5, 5, 4, 6, 2, 5, 6, 10, 9, 1, 6, 6, 4, 5, 4,
2, 4, 2])
hist * len(data) # Умножаем на количество чисел ввиду того, что мы подсчитывали долю, а не количество.
array([ 7., 7., 5., 5., 4., 6., 2., 5., 6., 10., 9., 1., 6.,
6., 4., 5., 4., 2., 4., 2.])
Визуальное сравнение
Можно сравнить визуально выведя обе гистограммы как график.
# Наивный способ отрисовки двух графиков столбцов.
plt.bar( elems, hist * len(data) ) # !len(data)
plt.bar( elems, hh[0] ) # Первый график был закрашен вторым.
<BarContainer object of 20 artists>
# Чтобы нарисовать два графика с колонками нужно изменить параметры отрисовки.
# Стандартно колонки имеют ширину 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(["Наш подход", "Батарейка"])
<matplotlib.legend.Legend at 0x7fe9dc6c44e0>
plt.bar( elems, hist * len(data), width = 0.4, align = 'center' )
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(["Наш подход", "Батарейка"])
<matplotlib.legend.Legend at 0x7fe9dc68a668>
Визаульно конечно можно получить некое общее представление, что все вопрядке, но хотелось бы более нажедного способа.
Можно по-элементно сравнить массивы на совпадение.
Численное сравнение
hh[0]
array([ 7, 7, 5, 5, 4, 6, 2, 5, 6, 10, 9, 1, 6, 6, 4, 5, 4,
2, 4, 2])
hist
array([0.07, 0.07, 0.05, 0.05, 0.04, 0.06, 0.02, 0.05, 0.06, 0.1 , 0.09,
0.01, 0.06, 0.06, 0.04, 0.05, 0.04, 0.02, 0.04, 0.02])
len(hist)
20
eq = hh[0] == hist*len(data)
eq # Опять же визуально видно, что все вроде истинно...
array([False, False, True, True, True, True, True, True, True,
True, True, True, True, True, True, True, True, True,
True, True])
np.array([True, True, True]).all(), np.array([True, False, True]).all()
(True, False)
np.array([False, False, False]).all()
False
eq.all( ) # Проверяет истинные ли все элементы вместе взятые.
False
eq == False
array([ True, True, False, False, False, False, False, False, False,
False, False, False, False, False, False, False, False, False,
False, False])
# Строим массив мест для которых ложь.
w = np.where( eq == False )
w
(array([0, 1]),)
hh[0][w], (hist*100)[w]
(array([7, 7]), array([7., 7.]))
(hist*100)[w]-hh[0][w]
array([8.8817842e-16, 8.8817842e-16])
np.allclose(hh[0], hist*len(data))
True
# Округляем числа до целого числа.
myhist = (hist*len(data)).round(decimals=14)
eq = hh[0] == myhist
eq
array([ True, True, True, True, True, True, True, True, True,
True, True, True, True, True, True, True, True, True,
True, True])
eq.all() # Теперь должна дать истину.
True
?np.round
Signature: np.round(a, decimals=0, out=None) Docstring: Round an array to the given number of decimals. See Also -------- around : equivalent function; see for details. File: /data/conda/anaconda3/envs/py36/lib/python3.6/site-packages/numpy/core/fromnumeric.py Type: function
(hh[0] - hist*len(data)).std() # Но можно и так, вычислив среднекадратичное отклонение разности.
2.664535259100376e-16
Упр. Как циклом за один найти наиболее часто встречаемый элемент?
Рассмотрим непрерывную случайную величину.
data = np.random.rand( 100 ) # Массив из 100 чисел равномерно распределеных на интервале [0 1).
np.min( data ), np.max( data ) # Убедимся в этом.
(0.029457044383584097, 0.9794707754246657)
В данном случае мы не сможем посчитать количества раз которое та или иная величина встречается. Для адекватного подсчета нужно разбить интервал на равные подинтервалы.
# Допустим мы хотим разбить интервал [0, 1) на десять равных подинтервалов.
d = np.linspace(0, 1, 11) # Тогда нам нужно на этом интервале расположить 11 точек.
d
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])
d[2], d[3] # Пример интервала.
(0.2, 0.30000000000000004)
get = d[2] <= data # Это даст те числа которые больше равны d[2], т.е. 0.2. gt -- greater or equal than.
lt = data < d[3] # А это те которые строго меньше 0.3. l -- less than.
get & lt # Объеденить их можно так.
array([False, False, False, True, False, False, False, False, False,
False, False, False, True, False, False, False, False, False,
False, True, 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, True, False, True, False, False, False, False,
False, False, False, False, False, False, False, False, True,
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, True, False,
False])
# Напишем функцию аналогичную дискретному случаю.
cnt = lambda data, e : np.mean( (e[0]<= data) & (data < e[1]) ) # Выборочная функция.
cnt( data, (d[2], d[3]) ), cnt( data, (d[5], d[6]) ) # Проведем подсчет для пары итервалов.
(0.1, 0.08)
interval = list(map( lambda i: (d[i], d[i+1]), range(10) ) ) # Создаем список кортежей границ интервалов.
interval # Содержит кортежи интервалов.
[(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)]
hist2 = map( lambda e: cnt( data, e), interval)
hist2 = list(hist2)
hist2
[0.05, 0.1, 0.1, 0.14, 0.05, 0.08, 0.14, 0.09, 0.13, 0.12]
plt.bar( d[:-1], hist2 ) # Какой то странный вид гистограммы для равномерного распределения....
<BarContainer object of 10 artists>
plt.bar( d[:-1]+0.05, hist2, width=0.09 ) # Изменена ширена колонки дабы они не накладывались друг на друга.
plt.xticks( d[:-1]+0.05 ); # 0.05 Дабы центр колонки был там где надо.
plt.hist( data ) # Опять же можно воспользоватся готовой функцией.
(array([ 8., 8., 10., 13., 6., 7., 12., 11., 10., 15.]),
array([0.02945704, 0.12445842, 0.21945979, 0.31446116, 0.40946254,
0.50446391, 0.59946528, 0.69446666, 0.78946803, 0.8844694 ,
0.97947078]),
<BarContainer object of 10 artists>)
Поиск объекта на изображении
img = plt.imread('7.1.02.tiff')
img.shape
(512, 512)
plt.imshow(img, cmap='gray')
<matplotlib.image.AxesImage at 0x7fe9dc212748>
sum0 = np.mean( img, axis = 0 )
sum1 = np.mean( img, axis = 1 )
sum0.shape, sum1.shape
((512,), (512,))
plt.plot( sum0 )
[<matplotlib.lines.Line2D at 0x7fe9dc17d240>]
plt.plot( sum1 )
[<matplotlib.lines.Line2D at 0x7fe9dc154d30>]
p = np.polyfit( np.arange(0, 512, 1), sum1, 1 )
p # p(x) = p[0] * x**deg + ... + p[deg]
array([4.66239766e-02, 1.63422512e+02])
xx = np.arange(0, 512, 1)
n = np.sum((xx - np.mean(xx))*(sum1 - np.mean( sum1 )))
d = np.sum((xx - np.mean(xx))**2)
n, d
(521478.361328125, 11184768.0)
plt.plot( sum1 - p[0]*xx)
[<matplotlib.lines.Line2D at 0x7fe9dc0bac88>]
np.min( sum0 )
155.58984375
xcrd = np.where( sum0 == np.min( sum0 ) )
xcrd
(array([244]),)
ycrd = np.where( sum1 == np.min( sum1 ) )
ycrd
(array([225]),)
plt.imshow(img, cmap='gray')
plt.plot( xcrd, ycrd, 'r*')
[<matplotlib.lines.Line2D at 0x7fe9dc0dba58>]
Возвращаемся к вопросу о формировании произвольного распределения по равномерному.
import scipy.stats as models # Загружаем модуль знающий о сложных распределениях.
import scipy.special as sp
x = np.random.rand(1000)
nn = sp.erfinv( 2*x - 1 ) * np.sqrt(2)
nn[20:40]
array([ 0.45951959, 0.15588191, -1.03389108, 0.91643848, -0.06141641,
-0.46014208, -0.38725496, 0.05105583, -0.38119611, 0.7527464 ,
1.11606978, -0.33193665, 0.67245078, 0.03220227, -0.50581301,
1.10607636, -1.76735953, 0.5705234 , 0.44561383, 0.34329685])
plt.hist( nn )
(array([ 1., 9., 42., 111., 226., 299., 207., 74., 23., 8.]),
array([-3.91441368, -3.18315763, -2.45190158, -1.72064553, -0.98938948,
-0.25813343, 0.47312262, 1.20437868, 1.93563473, 2.66689078,
3.39814683]),
<BarContainer object of 10 artists>)
# Можно приблизить плотность распределения за счет ядерного подхода.
gg = models.gaussian_kde( nn )
x = np.linspace(-3, 3, 100)
plt.plot( x, gg(x) )
[<matplotlib.lines.Line2D at 0x7fe9b6758390>]
1/np.sqrt(2*np.pi) # Максимальное значение плотности.
0.3989422804014327
np.max( gg(x) )
0.411074097771644
np.argmax( gg(x) )
49
gg(x[49])
array([0.4110741])
imax = np.argmax( gg(x) )
f'индекс {imax} в массиве x. В точке {x[imax]:.2} достигается максимальное значение {gg( x[imax] )[0]:2}'
'индекс 49 в массиве x. В точке -0.03 достигается максимальное значение 0.411074097771644'
gg = models.gaussian_kde( nn, 0.1 ) # Можно задать параметр
x = np.linspace(-3, 3, 100)
plt.plot( x, gg(x) )
[<matplotlib.lines.Line2D at 0x7fe9b66c96a0>]
gg = models.gaussian_kde( nn, 0.01 )
x = np.linspace(-3, 3, 100)
plt.plot( x, gg(x) )
[<matplotlib.lines.Line2D at 0x7fe9b66a5fd0>]