Введение в машинное обучение.
В дополнение к первой заметке рассматривются такие понятия как собственные функции, классы и иные библиотек относящихся к анализу данных.
Это предварительная версия! Любые замечания приветсвуются.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as models
import pandas as pd
# import seaborn
%config InlineBackend.figure_format = 'svg'
flat_df = pd.read_csv('berkeley_case.csv', sep = ';') # Считываем данные из файла.
flat_df # Выводим считанную таблицу. Она выглядит достаточно опрятно и эффектно.
flat_df.columns # Выводим название колонок.
# Переименовываем на русский лад.
flat_df.columns = ['Факультет', 'Пол', 'Параметр', 'Количество']
flat_df
# Выполним быструю агрегацию по полу.
total_stats = pd.pivot_table(flat_df, aggfunc = sum, index = 'Пол', columns = 'Параметр', values = 'Количество')
total_stats
total_stats.accepted/total_stats.applied # В процентах.
total_stats['Процент поступ'] = total_stats.accepted/total_stats.applied
total_stats
#total_stats['perc_admitted'] = [round_2digits( it ) for it in 100*total_stats.accepted/total_stats.applied ]
total_stats['Процент поступ'] = [ round( it, 2) for it in 100*total_stats.accepted/total_stats.applied ]
total_stats
df = pd.pivot_table(flat_df, index = 'Факультет', values = 'Количество', columns = ['Пол', 'Параметр'])
df
df['women']
df['women', 'applied']
df['accepted'] # так нельзя
idx = pd.IndexSlice
idx[:, 'accepted']
idx[:]
df.loc[idx[:]] # Никаких изменений.
df.loc[idx[:], idx[:, 'accepted']]
df_total = df['men'] + df['women']
df_total
df_totalT = df_total.T # Транпонируем таблицу.
df_totalT
df_totalT['пол'] = 'всего'
df_totalT
df_totalT.set_index('пол', append = True, inplace = True)
df_totalT
df_totalT = df_totalT.reorder_levels(['пол', 'Параметр'])
df_totalT
df_total = df_totalT.T
df_total
df = pd.concat([df, df_total], axis = 1)
df
df_inv = df.reorder_levels(['Параметр', 'Пол'] )
df_inv = df.reorder_levels(['Параметр', 'Пол'], axis = 1)
df_inv
df_inv = df_inv.sort_index(level = 0, axis = 1)
df_inv
100*df_inv.accepted/df_inv.applied # Вычисляем итоговые проценты на факультетах в зависимости от пола.
# Итоговые проценту получили, но их нужно округлить. Как раньше не получится. Объект сложный. Итерация не годится.
# Можно ввести округление.
# Для этого понядобятся функции
def round_2digits( x ): # Данная функция от одного аргумента. Как синус, квадратный корень и тому подобное.
return round(x, 2) # Вызываем функцию из самого Python, число 2 указывает на то что округляем две цифры после запятой.
round_2digits( 2.123456 )
round_2digits( 2.1 )
round_2digits( 2.12567 ) # Это именно округление, а не отрезание чисел после второго разряда после запятой.
admitted_perc = 100*df_inv.accepted/df_inv.applied
admitted_perc = round_2digits( admitted_perc )
admitted_perc
admitted_perc[['всего', 'men', 'women']].plot(kind = 'bar', title = 'Процент посупивших в Беркли')
df_applied = flat_df[ flat_df.Параметр == 'applied' ]
df_applied
gndr_fclt_appl = pd.pivot_table( df_applied, index = 'Факультет', values = 'Количество', columns = 'Пол')
gndr_fclt_appl
gndr_fclt_appl.sum()
gndr_fclt_appl = 100 * gndr_fclt_appl / gndr_fclt_appl.sum()
gndr_fclt_appl
round_2digits( gndr_fclt_appl )
faculty_stats = admitted_perc[['всего']].join( gndr_fclt_appl )
faculty_stats
faculty_stats.columns = ['Процент посупивших', 'доля мужчин', 'доля женщин']
faculty_stats.plot(kind = 'bar', title = 'Статистика факультетов Беркли')
Функции
def mysum( a, b): # Теперь у функции два аргумента
return a + b * 2; # Возвращаем значение выражения.
mysum( 3, 2) # 3 + 2 * 2
mysum( 3, a = 2) # При таком вызове система запутается. Точнее будет повторное присвоение первому аргументу.
mysum( b = 3, 2) # А так нельзя потомучто после присвоения переменной значения по имени позиции не учитываются.
mysum( b = 3, a = 2) # 2 + 3 * 2
# Можно сделать значения по умолчанию
def mysum2( a, b = 5):
return a + b * 2;
mysum2( 3 ) # 3 + 5 * 2
mysum2( 3, 2) # 3 + 2 * 2
# После переменных с значением по умолчанию обычные переменные не могут следовать
def mysum3( a, b = 5, c):
return a + b + c
# Тело функции ествественно может быть сложным
def doPow( a, c ): # Вычисляем возведение в степень. a возводится в степен c.
d = 1
for i in range(c):
d = d * a
return d
doPow( 2, 3)
doPow( 3, 4)
# Функция может иметь несколько точек воврата
def myAbs( a ):
if a < 0:
return -a # else опустили для краткости.
return a
myAbs( -11 )
doPow( 3, -5) # Для отрицательной степени не сработали.
# Можем вызвать сами себя:
def doPow( a, c ):
if c < 0: #Учитываем отрицательость степени
return doPow( 1/a, -c )
d = 1
for i in range(c):
d = d * a
return d
doPow( 3, -5)
(1/doPow( 3, -5))/3/3/3/3/3
# Можем и более сложные алгортмы делать. Например вычисление факториала.
def fact( n ):
if n <= 1:
return 1
return n * fact( n - 1 )
fact( 3 )
# Раз уж if был упомянут...
data = np.random.randint(-10, 10, 10) # Равномерное целочисленное распределение между -10 и 10 невключительно.
data
cnt_n = 0 # Количество отрицательных
cnt_p = 0 # Количество положительных
cnt_e = 0 # Количество равных нулю
data = np.random.randint(-10, 10, 100)
for d in data:
if d > 0:
cnt_p += 1 # cnt_p = cnt_p + 1
elif d < 0:
cnt_n += 1 # cnt_n = cnt_n + 1
else:
cnt_e += 1 # cnt_e = cnt_e + 1
print( ">0 : ", cnt_p, " <0 : ", cnt_n, " ==0 : ", cnt_e)
100/20 # Приблизительно равно количеству равных 0.
x = models.norm.rvs(2, 3, size = 100)
models.norm.fit( x )
np.mean( x )
np.std( x ) # Вот и все, т.е. по формулам.
pi = 0.3
lmbd = 3
def zip_pmf( x, pi = pi, lmbd = lmbd ):
if pi < 0 or pi > 1 or lmbd <= 0:
return np.zeros_like( x )
else:
return (x == 0) * pi + (1 - pi) * models.poisson.pmf( x, lmbd)
zip_pmf( 0 )
zip_pmf( 1 )
models.poisson.pmf( 0, lmbd)
fig, ax = plt.subplots(figsize=(8, 6))
xs = np.arange(0, 10)
palette = ['g','b']
ax.bar( 2.5*xs, models.poisson.pmf( xs, lmbd), width = 0.9, color = palette[0], label='Poisson')
ax.bar( 2.5*xs + 1, zip_pmf( xs ), width = 0.9, color = palette[1], label='Zero - Poisson')
ax.set_xticks( 2.5 * xs + 1 )
ax.set_xticklabels( xs )
ax.set_xlabel('$x$')
ax.set_ylabel('$P(X=x)$')
ax.legend()
N = 1000
iz = models.bernoulli.rvs( pi, size = N)
x = (1 - iz) * models.poisson.rvs( lmbd, size = N)
print(x[::50])
x.shape
fig, ax = plt.subplots(figsize=(8, 6))
ax.hist( x, width=0.8, bins = np.arange( x.max() + 1 ), normed = True)
ax.set_xticks( np.arange( x.max() + 1 ) + 0.4 )
ax.set_xticklabels( np.arange( x.max() + 1 ) )
ax.set_xlabel('$x$')
ax.set_ylabel('Частота')
from statsmodels.base.model import GenericLikelihoodModel
class ZeroInflated(GenericLikelihoodModel):
def __init__(self, endog, exog = None, **kwds):
if exog == None:
exog = np.zeros_like( endog )
super( ZeroInflated, self).__init__( endog, exog, **kwds)
def nloglikeobs( self, params):
pi = params[0]
lmbd = params[1]
return -np.log( zip_pmf( self.endog, pi = pi, lmbd = lmbd) )
def fit(self, start_params = None, maxiter = 10000, maxfun = 5000, **kwds):
if start_params is None:
lmbd_start = self.endog.mean()
extra_zeros = ( self.endog == 0 ).mean() - models.poisson.pmf(0, lmbd_start )
start_params = np.array([extra_zeros, lmbd_start])
return super( ZeroInflated, self ).fit(start_params = start_params, maxiter = maxiter, maxfun = maxfun, **kwds)
model = ZeroInflated( x )
res = model.fit()
res.params # Profit!
Используем распределение Бернули
Задается вероятностью p происхождения события. Например, вероятность кликнуть на банер, зайти на сайт, прийти в магазин или купить товар. Это все разные вероятности. Важное, что некое событие может произойти и неоднократно.
models.bernoulli.rvs( 0.7, 10) # Так нельзя, так как 10 будет использоватся как сдвиг.
В питоне функции могут использовать два типа аргументов или иначе параметров, которые ей передаются. Первый тип позиционый, который является стандартным для большинства языков программиования. Суть в том, что параметры присваиваются в порядке и задания при вызове. А бывают аргументы именные, которые можно присвоить вне очереди.
Рассмотреть документацию по данной функции https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.bernoulli.html
Конкретно изучим описание функции rvs: rvs(p, loc=0, size=1, random_state=None)
Последнее к сожалению позволит понять какие значения были присвоены аргументам по умолчанию. часть аргументов просто перечислена (а данном случае только p), а часть через знак '='.
# Вообще у любого объекта, функции по соглашению может быть поле __doc__ которое дает краткую документацию по нему.
# Так как это поле является строчкой, то для аккуратного вывода лучше применить функцию print
# Подробнее о ней будет позже. Сейчас она используется для вывода самоописания функции.
print( models.bernoulli.rvs.__doc__ )
Из данного описания следует, что сначала идут позиционные параметры (arg1, arg2, ...), а потом именные (loc, size, random_state). Тем не менее, если они идут подряд, то присвоение идет подряд.
# В данном случае:
models.bernoulli.rvs( 0.7, 10) # 0.7 это параметр распределения, а 10 присвоилось аргументу loc.
models.bernoulli.rvs( 0.7, 10, 5) # Вот теперь очередь дошла до размера, которому было присвоено значение 5.
# Но можно было бы и явно присвоить переменной отвечащей за размер размер:
models.bernoulli.rvs( 0.7, size = 10 ) # Аргументу size присвоено значение вне очереди.
# Нельзя позиционный аргумент подавать после именных. Последнее выдаст соответствующую ошибку:
models.bernoulli.rvs( size = 10, 0.7 )
# Но если аргументы заданы по имени, то их можно подавать функции в произвольном порядке.
models.bernoulli.rvs( size = 10, p = 0.7 )
models.bernoulli.rvs( size = 10, p = 0.7, loc = 5 ) # В данном случае было задано и значение параметра loc, положения.
Закон больших чисел
Нам нужно взять среднее количество произошедших событий. Точнее понять как это среднее рапределено. Теория говрит, что оно должно стремится к нормальному распределению.
# Имитация закона больших чисел:
data = [] # Содали пустой список
for i in range( 1000 ): # Мы явно указали размер нешей выборки
lap = models.bernoulli.rvs( 0.7, size=100 ) # а также количество элементов в устреднение.
data.append( np.sum( lap )/100 )
plt.hist( data );
# Хотим цикл ещё раз запустить с другим количеством эллементов в сумме.
Продолжаем с Бернули
# Вводим понятие функции, т.е. сами пишем.
def buildBernSum( p, cycles = 1000, elems = 100 ): # Два аргумента. Оба необязательны. Имеют значение по умолчанию.
data = [] # Содали пустой список
for i in range( cycles ):
lap = models.bernoulli.rvs( p, size=elems ) # 0.7
data.append( np.sum( lap )/elems )
return data # Вернули результат.
bern = buildBernSum( 0.7 )
plt.hist( bern );
bern = buildBernSum( 0.7, elems = 10000 )
plt.hist( bern );
bern = buildBernSum( elems = 10000, p = 0.7)
plt.hist( bern );
bern = buildBernSum( 0.7, 100000, 10000 ) # Позиционное присвоение параметров
plt.hist( bern );
bern = buildBernSum( 0.7, 4000, 3000 )
plt.hist( bern );
Статистика Шапиро-Уилка
np.random.seed( 2018 ) # Инициализация генератора случайных чисел.
dat = models.norm.rvs( 2, 3, 10 ) # 2 и 3 это параметры нормального распределения.
plt.hist( dat );
t = models.shapiro( dat ) # Шапиро-Уилка тест для нормальности. Предполагается, что данных не очень много ( < 5000).
t
# t[0] -- Это значение самой статистики.
# t[1] -- Это значение p-value
print( "Статистика -- ", t[0], ", p-значение ", t[1])
Смотрим на p-значение. Если оно больше, например, 5%, то говорим, что с такой значимостью не можем отвергнуть гипотезу о нормальности даных.
plt.hist( dat * dat ); # Строим гистограмму для квадратов нормального распределения. Для данного распределения
models.shapiro( dat * dat ) # p-значение очень маленькое. Поэтому гипотезна о нормальности отвергнута.
plt.hist( 1/dat ); # Строим гистограмму для обратной величины от нормального распределения.
models.shapiro( 1/dat ) # Она тем более не являются нормальными. p-значение ну очень маленькое.
dat = models.norm.rvs( 2, 3, 10 ) # Возьмем новый набор данных.
plt.hist( dat * dat );
models.shapiro( dat * dat ) # p-значение чуть повыше, но все-равно меньше традиционной значимости 0.05.
plt.hist( 1/dat ); # При молой выборке конечно возможны ложные результаты
models.shapiro( 1/dat ) # Действительно смахивает на нормальное распределение. p-значение значительно выше 5%.
dat = models.norm.rvs( 2, 3, 100 ) # Если взять больше данных, то вероятность такого мала.
plt.hist( dat ); # С увеличением данных, статестический тест может давать странные результаты.
models.shapiro( dat ) # В данном случае хоть тест и пройден, все-таки 8% не сильно далеко от 5%.
plt.hist( dat * dat ); # Для произведения нормальных распределений.
models.shapiro( dat * dat ) # p-значение совсем маленькое. Поэтому, тест не пройден.
plt.hist( 1/dat ); # При большой выборке вероятность ложного результата для сильных отклонений от нормального мала.
models.shapiro( 1/dat ) # p-значение совсем маленькое. Поэтому, тест не пройден.
dat = models.uniform.rvs( size = 10 ) # Возьмем проcто-напросто другой тип раcпределения, а именно равномерное.
plt.hist( dat );
models.shapiro( dat ) # Тем не менее тест ломается на равномерном распределении, так как p-значение пройдено.
dat = models.uniform.rvs( size = 100 ) # Теперь увеличим объем выборки для равномерного распределения.
plt.hist( dat );
models.shapiro( dat ) # p-значение не пройдено, поэтому гипотезу отвергаем.
dat = models.f.rvs( 5, 7, size = 100 ) # f распределение с параметрами 5 и 7.
plt.hist( dat );
models.shapiro( dat ) # f распределении отвергаем, так как p-значение существенно меньше.
dat = models.f.rvs( 20, 35, size = 100 ) # f распределение с параметрами 20 и 35.
plt.hist( dat );
models.shapiro( dat ) # Для данных параметров оно тоже отвергнуто.
dat = models.norm.rvs( 2, 3, 1000 )
plt.hist( dat );
models.shapiro( dat ) # Все хорошо с p-значением.
#С увеличением параметров распределение будет стремится к нормальному
Anderson-Darling
Предыдущая статистика работает только для нормального распределения. Для данного же теста можно задать используемое распределение: norm, expon, logistic, gumbel, gumbel_l, gumbel_r, extreme1. Пока будем оставатся в рамках нормального.
В Шапиро-Уилка граница была единственной. Той что обычно принято пользоваться. В методе же Андерсона таких границ считается несколько для различных уровней значимости.
dat = models.uniform.rvs( size = 100 ) # Пвторно проверим на равномером распределении.
plt.hist( dat );
models.anderson( dat, dist='norm' ) # Тест отвергает данные как принадлежащие равномерному распределению.
# Статистику нужно возвести в квадрат и сравнить с критерием. Если больше, то гипотеза отвергается.
Проверка того что распределение то что нужно.
dat = models.norm.rvs( 2, 3, 100 ) # Для начала прверим для нормального распределения.
res = models.anderson( dat )
res
res[1] # Значение значимости критических уровней.
res[2] # Сами уровни.
# Можно представить их матрицей.
table = np.array([res[1], res[2]])
table
# Хотим красивое представление результата/критерия в виде таблицы
import pandas as pd # Очередная библиотека полезная для представления данных в таблицах, и их считывания.
pd.DataFrame( table, columns=['Уровни', 'Значимость'] ) # Система ругнется так как матрица не правильного размера.
# Данная функция ожидала матрицу состоящую из двух столбцов.
# В Numpy естественно есть функция для транспонирования матрицы.
np.transpose( table )
# Транспонироваие дает результат и в общем случае: она переставляет оси.
x = np.ones((1, 2, 3)) # ones вспомогательная функция составляющая многомерную матрицу из единиц. В скобках указан размер.
x
xT = np.transpose(x, (1, 2, 0)) # В скобках указана перестановка исходных осей, т.е. оси будут идти в порядке 1, 2, 0.
xT
x.shape, xT.shape # Если оси переставлены, то мы это увидем на соответствующем образом изменившимся порядке их размере.
# Вначале должна идти первая компонента, потом вторая, а уже потом нулевая векто рараеров (1, 2, 3).
# Тогда получим 2, 3, 1
# Тогда:
pd.DataFrame( np.transpose( table ), columns=['Значение статистики', 'Уровень значимости'] )
# Получаем красивую таблицу.
# Но можно иначе.
table = np.dstack( (res[1], res[2]) )
table, table.shape # Но появится лишняя размерность.
# Её можно убрать взяв первый элемент:
table[0], table[0].shape
# тогда опять можем получить таблицу.
pd.DataFrame( table[0], columns=['Значение статистики', 'Уровень значимости'] )
# Скорее всего ошибка в коде. Значимости нужны в обратном порядке.
table = np.dstack( (res[1], res[2][::-1] ) )[0]
pd.DataFrame( table, columns=['Значение статистики', 'Уровень значимости'] )
# В res[0] у нас значение статистики.
res[0] # Для проверки гипотезы нужно сравнить это число с критическими значениями.
Бернули
bern = buildBernSum( 0.7, 4000, 3000 )
res = models.anderson( bern ) # Второй параметр указывает на распределение. Поумолчанию проверка на нормальность.
res # Дествительно стремимся к нормальному распределению.
params = models.norm.fit( bern )
params
models.norm.interval( 0.95, params[0], params[1] )
models.bernoulli.interval( 0.7, 0.7 ) # Системе плохо.
models.norm.interval( 0.95 )
pp = np.mean( bern )
pp
dd = np.sqrt(pp*(1-pp)/3000)
dd
lim = models.norm.interval( 0.95, 0, 1 )
lim
pp + dd*lim[0], pp + dd*lim[1]
Pearson
models.pearsonr([1, 2, 3, 4, 5], [4, 7, 10, 13, 16])
plt.plot( [1, 2, 3, 4, 5] )
plt.plot( [4, 7, 10, 13, 16] )
models.pearsonr([1, 2, 3, 4, 5], [-12, -14, -16, -18, -20])
models.pearsonr([1, 2, 1, 3, 1], [10, 12, 10, 14, 10]) # Необязательно монотонный данные.
plt.plot( [1, 2, 1, 3, 1] )
plt.plot( [10, 12, 10, 14, 10] )
Spearman
models.spearmanr([1,2,3,4,5], [50,51,60,100,200])
plt.plot( [1,2,3,4,5] )
plt.plot( [50,51,60,100,200] )
models.spearmanr([1,2,6,4,5], [50,51,300,100,200])
plt.plot( [1,2,6,4,5] )
plt.plot( [50,51,300,100,200] )
models.spearmanr([1,5,4,3,2], [5,6,7,8,7])
models.laplace.rvs()
A/B Тест