Механико-математический факультет МГУ
Магистратура по направлению "Искусственный интеллект" 2022-23

Весенний семестр

10 Feb 2023

Разработка оконных приложений на Питоне
с помощью модуля tkinter

Задача "построение графика интерполяционного многочлена
по заданным узлам интерполяции"

Задача линейной регрессии и ее решение методом
наименьших квадратов

Задача интерполяции

Даны узлы интерполяции
    x0, x1, x2, ..., x_n       x_i != x_j
    y0, y1, y2, ..., y_n
Существует единственный многочлен степени n
    f(x) = a0*x^n + a1*x^(n-1) + ... + a_n
такой, что
    f(x_i) = y_i, i = 0, 1, ..., n
(имеем систему с матрицей Вандермонда, определитель которой
отличен от нуля)
    x0^n x0^(n-1) ... 1
    x1^n x1^(n-1) ... 1
    . . . 
    xn^n xn^(n-1) ... 1

Многочлен в форме Ньютона
    p_n(x) = a0 + a1(x-x0) + a2(x-x0)(x-x1) + ... + a_n(x-x0)...(x-x_n-1)  
Добавим узел x_n+1, y_n+1
    p_n+1(x) = p_n(x) + a_n+1(x-x0)(x-x1)...(x-x_n-1)(x-x_n)
    a_n+1 = (p_n+1(x) - p_n(x))/( (x-x0)(x-x1)...(x-x_n-1)(x-x_n) )
подставим x = x+1
    a_n+1 = (y_n+1 - p_n(x_n+1)/( (x_n+1-x0)(x_n+1-x1)...(x_n+1-x_n-1)(x_n+1-x_n) )
    
        ПЕРЕРЫВ 10 минут до 18:30
        
Задача линейной регрессии

Есть две величины X, Y, X -- независима переменная, Y -- зависимая
Хотим построить модель зависимости между X и Y

Модель определяется набором параметров w0, w1, ..., w_n
Зависимость от параметров линейная
    Y = w0*f0(X) + w1*f1(X) + ... + w_n*f_n(X)
В частности, рассмотрим задачу нахождения наилучшего интерполяционного
многочлена
    Y = w0*x^n + w1*x^n-1 + ... + w_n = p(x)
пусть заранее известно, что степень многочлена deg(p) = n
и нам известны значения Y в m точках, где m >= n

Если m > n, то точного решения в общем случае нет,
попробуем найти решение, минимизирующее ошибку.
В качестве функции ошибки берется сумма квадратов разностей
(квадрат евклидова расстояния между [p(X_i)] и [Y_i]
для тренировочных данных [X_i, Y_i]
    
    w0*X_0^n + w1*X_0^n-1 + ... + w_n = Y_0        
    w0*X_1^n + w1*X_1^n-1 + ... + w_n = Y_1        
    . . .
    w0*X_m^n + w1*X_m^n-1 + ... + w_n = Y_m       
    
Несовместная система для m > n

    AW = Y
    A: R^n --> R^m
    AW -- это линейное подпространство в R^m
   

A -- матрица Вандермонда
    X0^n  X0^n-1 ... 1
    X1^n  X1^n-1 ... 1
    . . .
    X_m^n X_m^n-1 ... 1
    
Попробуем найти такое решение  W = V, что квадрат расстояния
между AV и Y миимален.

Понятно, что V -- это основание перпендикуляра,
опущенного из Y на подпространство AW, 
то есть вектор AV - Y перпендикулярен подпространству AW

    (AV - Y, AW) = 0
    (AV - Y)^T * AW = 0
 Так как это выполнено для любого w <- R^n, то
    (AV - Y)^T * A = 0
    V'A'A - Y'A = 0 (через штрих обозначено транспонирование)
    V'A'A = Y'A   транспонируем обе части равенства
    A' A V = A' Y
Это нормальное уравнение в методе наименьших квадратов
У матрицы A столбцы линейно независимы =>
матрица A'A обратима =>
    V = (A' A)^(-1) A' Y
    V = inverse(A' A) A' Y
Определение:
матрица
    pinv(A) = inverse(A' A) A'
называется псевдообратной к матрице A
Решение задачи в методе наименьших квадратов:
    V = pinv(A)*Y

===============================================

10 Mar 2023

    Метод опорных векторов
    
Мы хотим построить линейный классификатор

Наши объекты задаются набором параметров
       x = (x_1, x_2, ..., x_n) <- R^n
Линейный классификатор:
       По x мы вычисляем вещественное значение:
       <w, x> = w1*x1 + w2*x2 + ... + w_n*x_n
       w -- вектор весов
    f(x) = <w, x> - b
Если f(x) > 0, то мы относим объект x к классу 1,
иначе относим объект к классу 2.

Нужно построить классификатор f(x), то есть найти w <- R^n,
b <- R.

Имеется тренировочная база
     (x1, y1), (x2, y2), ..., (x_m, y_m),
     x_i <- R^n,   y_i = +-1
     
     y_i = +1, если x_i принадлежит классу 1,
     y_i = -1, если x_i принадлежит классу 2.
     
Как работает классификатор?
     f(x) > 0  ==>  относим x к классу 1,
     f(x) < 0  ==>  относим x к классу 2.
     
Гиперплоскость f(x) = 0 должна максимально надежно
разделять два класса точек (из тренировочной базы).
В разделимом случае это означает, что минимальное расстояние
до точек первого класса равно мин. расстоянию до точек второго
класса.

Предложение
Расстояние от точки x до гиперплоскости, заданной уравнением
    <w, x> - b = 0,
равно
    h = (<w, x> - b)/|w|
    
Доказательство
    Пусть w' = w/|w|
    x0 -- любая точка на гиперплоскости
    h -- расстояние от x до гиперплоскости
    h = <(x - x0), w'> = <x, w'> - <x0, w'> =
      = <x, w>/|w| - <x0, w>/|w| = <x, w>/|w| - b/|w| =
      = (<w, x> - b)/|w|
--      

Разделяющая гиперплоскость зависит только от опорных точек,
т.е. точек, модуль расстояния от которых
до гиперплоскости в точности равен h
Точка x_i опорная:
    y_i*(<w, x_i> - b)/|w| = h
Для остальных точек      
    y_i*(<w, x_i> - b)/|w| > h
В любом случае выполняется неравенство
    y_i*(<w, x_i> - b)/|w| >= h
Домножим это неравенство на 1/h, чтобы в правой части получить 1:
    y_i*(<w, x_i> - b)/(|w|*h) >= 1

Обозначим w' = w/(|w|*h), b' = b/(|w|*h)
Получим неравенство
    y_i*(<w', x_i> - b') >= 1
Так как w' = w/(|w|*h), то
    |w'| = |w|/(|w|*h) = 1/h
    h = 1/|w'|
Опустим штрих, получаем задачу:    
    y_i*(<w, x_i> - b) >= 1
    h --> max
Это эквивалентно
    |w| --> min
Это всё равно, что минимизировать квадрат нормы:
    y_i*(<w, x_i> - b) >= 1  для i = 1, 2, ..., m
    |w|^2 --> min
Решая эту задачу, мы найдем гиперплоскость (линейный
классификатор), разделяющую два класса точек.        
        
    НЕРАЗДЕЛИМЫЙ СЛУЧАЙ
    
В общем случае система неравенств           
    y_i*(<w, x_i> - b) >= 1  для i = 1, 2, ..., m
решения не имеет. Ослабим эти условия, позволим
неравенству быть слабее на величину ksi_i >= 0:
    y_i*(<w, x_i> - b) >= 1 - ksi_i,
                 ksi_i >= 0
ksi_i -- величина штрафа за то, что неравенство не
выполняется точно. Мы хотим найти минимальный штраф,
при котором система всё-таки имеет решение.

Мы получаем следующую задачу:    
    y_i*(<w, x_i> - b) >= 1 - ksi_i,    i = 1, 2, ..., m
    ksi_i >= 0
При этом мы хотим минимизировать ширину полосы,
разделяющей два класса точек, плюс суммарный штраф:
    |w|^2  + С * Sum_i ksi_i --> min

Давайте избавимся от ослабляющих переменных ksi_i.        
    ksi_i >= 1 - y_i*(<w, x_i> - b)
    ksi_i >= 0
Надо, чтобы ksi_i были как можно меньше ==>
    ksi_i = max(0, 1 - y_i*(<w, x_i> - b))

Мы свели нашу задачу к следующей задаче безусловной
минимизации: 
    |w|^2  + С * Sum_i max(0, 1 - y_i*(<w, x_i> - b)) --> min
Обозначим z = 1 - y_i*(<w, x_i> - b)
Функция
        H(z) = max(0, 1 - z)
называется hinge loss. Мы получаем следующую задачу:
    |w|^2  + С * Sum_i H(y_i*(<w, x_i> - b)) --> min
Нужно найти минимум этого выражения по переменным w, b.
        
=================================================

17 Mar 2023

Имеем 2 класса объектов, объект задается набором 
из n чисел, т.е. точкой (вектором) n-мерного пространства 
R^n. Нам нужно построить линейный классификатор:
    f(x) = <w, x> - b
    w -- вектор весов, w <- R^n
    x -- объект, x <- R^n,
    c -- свободный член, c <- R
    
Тренировочная база:
    x1, x2, ..., x_m
Разметка
    y1, y2, ..., y_m    
    y_i = { 1,  если x_i принадлежит классу 1,
          { -1, если x_i принадлежит классу 2.
          
Разделющая гиперплоскость задается уравнением
   f(x) = <w, x> - b = 0
Классификатор:
   если f(x) > 0, то объект x относим к первому классу,
                  иначе ко второму.ы   
              
==============================================

24 Mar 2023

Ядро K(x1, x2)

    phi: R^n --> R^k,   k > n или k = infinity
Отображение phi называется спрямляющим отображением
    K(x1, x2) = <phi(x1), phi(x2)>
    
Важно, что знать отображение phi и вычислять его не обязательно!
K(x1, x2) может просто задаваться в виде некоторой функции.
Не всякая функция является ядром, но есть набор действий,
которые строят новые ядра из уже известных. Нам важно,
что следующие функции явлются ядрами:
1) Гауссовское ядро
    rbf(x1, x2) = exp(-||x1 - x2||^2 / 2*sigma^2)                  
    Radial Basis Function
2) полиномиальное ядро степени d
    poly(x1, x2) = (<x1, x2> + c0)^d
3) сигмоидное ядро
    . . .
    
Домашнее задание:
применить ядровый метод опорных векторов
к задаче классификации рукописных цифр.
Определить наилучшие параметры среди
1) ядро -- linear, poly или rbf;
2) константа регуляризации C:
   [0.25, 0.5, 1., 2., 4.]
3) для ядра poly -- степень (degree):
   [2, 3, 4]
   
Цифры можно загрузить из модуля
from sklearn.datasets import load_digits

Цифры -- это матрицы 8x8, содержащие числа от 0 до 4,
то есть наши объекты имеют размерность 64.

Разметка объектов target[] -- числа от 0 до 9 (номера цифр).

В методе опорных векторов мы разделяем точки только двух классов
   y[i] = +-1
   
Когда у нас много классов, то для классификации можно
использовать либо метод "один против остальных"
    ovr  -- one versus rest,
либо "каждый против каждого"
    ovo  -- one versus one.     
В методе "ovo" как бы происходит круговой турнир.
Если число классов n, то строятся 
n*(n-1)/2 классификаторов, при предсказании
класса у конкретной точки мы применяем все классификаторы
clf[a, b] для пар классов a < b. Если классификатор
clf[a, b](x) при фиксированном классе a для всех b выигрывает 
максимальное число сравнений по всем 
классам a, то точка x относится к классу a.
(Аналогия с круговым турниром баскетбольных команд.)
    predict(x) = argmax_a { Sum_b clf[a, b](x) }
    
===================================================
31 Mar 2023


Борьба с выбросами в задаче линейной регрессии

    x1, x2, ..., x_m <- R^n -- значения независимой переменной
Имеется линейная функция
    y = f(x) = w1*x1 + w2*x2 + ... + w_n*x_n + b
    f(x) = <w, x> + b
Можно опустить b, если мы заменим вектор x <- R^n на
вектор x' = (x, 1) <- R^n+1
       w' = (w, b)
       f(x) = <w, x> + b = <w', b'>
Свободный член не используется.

Итак, рассматриваем задачу нахождения линейной зависимости
между независимой переменной x <- R^n и зависимой переменной
y <- R.
    y = <w, x>
Нужно найти вектор w.
На тренировочном наборе данных x1, x2, ..., x_m <- R^n
функция принимает известные значения y1, y2, ..., y_m.
В методе наименьших квадратов мы минимизируем функцию ошибки
(функцию потерь)
    L(w) = Sum ||<w, x_i> - y_i||^2     "MSE"
Проблема в том, что на эту функцию очень сильно влияют случайные
большие отклонения ("выбросы", outlayers).

В методе наим. квадратов мы минимизируем квадратичные отклонения

    L(w) = Sum E(<w, x_i> - y_i)     "MSE"
           E(z) = ||z||^2  -- выбросы сильно влияют
Хотелось бы, чтобы 
           E(z) = ||z||    -- линейная функция, выбросы
                              на результат почти не влияют              

Функция Хубера H(z) обладает свойствами этих обеих
функций: вблизи нуля она ведет себя как квадрат,
         для больших значений z -- как модуль z (линейная ф-ция). 
H(z) = { 1/2*z^2,   |z| <= 1
       { |z| - 1/2, |z| > 1


Таким образом, в задаче линейной регрессии можно минимизировать
функцию потерь на основе функции Хубера:
    L(w) = Sum H(<w, x_i> - y_i)     "Huber"
В методе наим. квадратов мы находим точное решение с
помощью псевдообратной матрицы. В методе Хубера
точного решения не существует, но можно минимизировать
функцию потерь разными стандартными методами минимизации
(метод сопряженных градиентов и т.п.).

В модуле sklearn.linear_model
имеются классы LinearRegression (для метода наим. квадратов
MSE) и HuberRegressor (при использовании функции Хубера).

Можно ли минимизировать сумму не квадратов, а просто модулей
отклонений? MAE -- Maximal Absolute Error
    L(w) = Sum |<w, x_i> - y_i|
    
Это можно сделать, используя линейное программирование.
В n-мерном пространстве задано выпуклое множество,
равное пересечению конечного количества полупространств.
Также задана целевая линейная функция. Нужно найти точку,
в которой эта функция принимает минимальное/максимальное
значение.

Проблема в том, что в нашей функции потери используется модуль.
Как быть?
Есть трюк, который позволяет избавиться от модуля.
    f = u - v,     u >= 0, v >= 0
    |f| <= u + v
При этом, если мы минимизируем сумму u + v, то |f| = u + v
    f = u - v,     u >= 0, v >= 0
        u + v --> min
=> либо u = 0, либо v = 0 => |f| = u + v

Вернемся к нашей задаче:
    L(w) = Sum |<w, x_i> - y_i| --> min
Введем переменные u_i, v_i и соотношения
    <w, x_i> - y_i = u_i - v_i,   i = 1, 2, ..., m
    u_i >= 0,  v_i >= 0
Целевая функция, которую мы минимизируем, равна
    Sum_i (u_i + v_i)  --> min
    
Итак, мы получили задачу линейного программирования:    
    <w, x_i> - y_i = u_i - v_i,   i = 1, 2, ..., m
    u_i >= 0,  v_i >= 0
    Sum_i (u_i + v_i)  --> min

В Питоне есть модуль pulp для задач линейного программирования
В задаче (объект myTask)
линейного программирования используются:
1) переменные
    -- свободные
        x = LpVariable("x")
    -- ограниченные снизу
        x = LpVariable("x", 0.)
    -- ограниченные сверху
        x = LpVariable("x", upBound = 10.)        
    -- ограниченные с обеих сторон
        x = LpVariable("x", lowBound = -1., upBound = 10.)
2) можно создать массив переменных
        w = LpVariable.dicts("w", range(n), lowBound = 0)
    создаются переменные w[i], i = 0,1,..., n-1
    с именами w_0, w_1, ..., w_n-1
3) можно указать ограничения
    myTask += 2*x + 3 <= 0
    myTask += x - 2*y == 1
4) массивы переменных и выражений можно суммировать
    myTask += lpSum([w[i] for i in range(n)]) > 1
5) целевая функция ровно одна (нет операции сравнения)
    myTask += x + 3*y
    или
    myTask += lpSum([u[i] + v[i] for i in range(n)])
6) для решения задачи вызываем метод solve()
    myTask.solve()
    статус решения:
    myTask.status
    
    
Решением являются значения переменных
    value(x)
    value(y)
    value(w[i])
    
==========================================

7 Apr 2023
Метод главных компонент
Principal Component Analysis (PCA)

==========================================

14 Apr 2023

Байесовский классификатор

Обучение с учителем

Дан набор объектов n-мерного пространства
и разметка этих объектов: для каждого объекта известно
значение функции или класса, к которому этот объект
принадлежит. Надо построить математическую модель,
которая по произвольному объекту вычисляет функцию / 
класс, к которому он должен принадлежать.

Байесовский классификатор исходит из того, что для каждого
класса нам известна плотность вероятности распределения
объектов этого класса в пространстве.

    P_i(x) = плотность вероятности распределения 
             объектов i-го класса
             
Пусть есть k классов и априори известны законы распределения
для всех этих классов. Тогда байесовский классификатор
выдает предсказание:
    predict(x) = argmax_i=1..k( P_i(x) )
    
Рассмотрим 3 варианта байесовских классификаторов:
1) наивный байесовский классификатор --
   нормальное распределение, причем координаты не зависят
   друг от друга (матрица ковариаций диагональная);
2) линейный классификатор LinearDiscriminantAnalysis
   каждого класса нормальное распределение, причем матрицы
   ковариаций для всех классов одинаковы;
3) квадратичный классификатор QadraticDiscriminantAnalysis   
   наиболее общий. Единственное ограничение -- 
   точки одного класса распределены по нормальному закону,
   параметры распределений могут быть разными для разных
   классов.


==========================================

21 Apr 2023
    
Кластеризация (обучение без учителя)

Есть набор объектов без разметки.
Мы хотим разбить их на классы, а также построить
классификатор, который для любой точки предсказывает
номер класса, к которому ее следует отнести
(метод predict).

Обучение классификатора осуществляется с помощью метода
fit.

Мы рассмотрим самые простые и популярные методы
кластеризации KMeans и DBSCAN.

KMeans основан на нахождении k средних значений (mean --
среднее значение в математике).

DBSCAN основан на нахождении областей с максимальной
плотностью объектов (рассмотрим следующий раз).

В методе нам априори известно число кластеров,
на которые мы должны разбить всё множество объектов
    x1, x2, ..., x_m,   x_i <- R^n, i = 1..m
KMeans основан на минимизации суммы дисперсий объектов
в каждом кластере C_i:
выбираем центры кластеров (центроиды) s1, s2, ..., s_k
так, чтобы сумма дисперсий объектов по всем классам
была бы минимальной:
    V_i = |C_i| * Sum       ||x_j - s_i||^2
                  x_j <- C_i
                  
Центроиды s1, s2, ..., s_k выбираются так, чтобы сумма
дисперсий была бы минимальной:

    (V1 + V2 + ... + V_k) --> min  по  s1, s2, ..., s_k
                  
Если центры кластеров (точки s1, s2, ..., s_k) фиксированы,
то любая точка x относится к тому кластеру, расстояние до
центра которого минимально.

Самый простой алгоритм кластеризации типа KMeans
следующий:
1) выбираем случайно k точек (например, любые случайные
   k точек среди объектов x1, x2, ..., x_m),
   и считаем их центрами кластеров s1, s2, ..., s_k;
2) разбиваем все точки на k кластеров, относя точку
   к тому кластеру, расстояние до центра которого минимально;
3) находим новые центроиды кластеров s1', s2', ..., s_k' и
   переходим обратно к пункту 2.
Повторяем эти две процедуры, пока центроиды s1, ..., s_k
меняются. 
Алгоритм заканчивается, когда на очередном шаге
    ||s_i - s_i'|| < eps  для всех i.

родемонстрируем этот алгорим на Python3.

import sklearn.cluster import KMeans

clf = KMeans(n_cluster = 3, n_init = 10)

clf.fit(x)  # Обучаем классификатор на объектах x (массив)

y_predicted = clf.predict(x)

======================================

28 Apr 2023

Методы кластеризации

Метод кластеризации K-Means предполагает, что априори
известно число классов. Дается множество точек n-мерного
пространства и число K классов, на которые эти точки надо
разбить. Результатом работы этого метода будет вычисленный набор
из K центроидов классов. Точка относится к тому классу,
центроид которого расположен ближе всех.

Недостатки:
1) вычисляется локальный минимум функции ошибки,
   при разных запусках могут получаться разные результаты;
2) хорошо работает только тогда, когда выпуклые оолочки
   разных кластеров не пересекаются.
   
Более интересный метод (самый популярный метод кластеризации)
это DBSCAN (Density-Based Spatial Clustering 
            of Applications with Noise)
Метод основан на пространственном анализе плотности точек,
при этом мы считаем, что могут быть выбросы, неправильно
измеренные объекты, шум и т.п., поэтому мы классифицируем
не все точки. Часть точек мы отмечаем как шум (noise) и не относим
их ни к какому кластеру.

Метод основан на двух переметрах:
    eps -- задает размер окрестности (по умолчанию
           используется евклидова метрика);
    minPts -- минимальное количество точек в eps-окрестности.
    
Определение: точка p считается основной точкой
             (core point), если в ее eps-окрестности
             имеется не меньше minPts точек, включая и
             ее саму.
             
Точка t считаеся достижимой (reachable) из основной точки p,
если t входит в eps-окрестность точки p.

Две точки t1 и t2 связные, если их можно связать цепочкой
основных точек:
    p1, p2, ..., p_k -- основные точки, причем
    t1 <- U(p1), где U(p1) -- eps-окрестность точки p1,
    t2 <- U(p_k),
    p2 <- U(p1), p3 <- U(p2),  ..., p_k <- U(p_k-1).
         
Кластер состоит из всех попарно связных точек.

Таким образом, мы разбиваем все исходное множество точек на
кластеры, а также на выбросы (шум) -- точки, которые 
не попали ни в один из кластеров.

В этом методе у кластеризатора нет метода predict.

from sklearn.cluster import KMeans, DBSCAN

clf_dbscan = DBSCAN(eps = 0.5, min_samples = 4)
clf_dbscan.fit(x)   # Выполняется алгоритм кластеризации

После этого имеем:

clf_dbscan.label_   массив номеров кластеров
                    Для i-й точки из массива x это будет
                    номер кластера, к которому эта точка
                    отнесена. Если точка не отнесена ни к
                    какому кластеру (шум, выброс), то
                    clf_dbscan.label_[i] = -1
                    
Можно также узнать, какие точки являются основными (core points):

clf_dbscan.core_sample_indices_  массив индексов основых точек.                    

============================================================

5 May 2023

Решающие деревья (Decision Trees)

С помощью решающих деревьев можно решать две задачи:
   -- задача классификации;
   -- задача регрессии
(обучение с учителем).

Дана база объектов X = x1, x2, ..., x_m
Разметка           y = y1, y2, ..., y_m.

В задаче классификации y_i -- номер класса, к которому
                              принадлежит объект x_i
                              
В задаче регрессии y_i <- R (вещественные числа).

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

Идея использования дерева:
    используем предикат, который разбивает множество
    объектов на 2 непересекающихся подмножества,
    причем разнообразие объектов внутри этих подмножеств
    уменьшается.
    
Пример экзамена по теории чисел в применении к криптографии:

                Малая теорема Ферма?
                 нет   /                \  да
               
              Что такое Z_m             вероятностный тест 
                                             простоты
               нет /      \ да                  /    \ да
                                               /    оценка 5
             оценка 2   вычислить обратный  оценка 4
                        элемент в Z_m
                         нет /      \ да
                          оценка 2  оценка 3
                          
Решающее дерево бинарное.

Терминальным вершинам (листьям)
соответствует номер класса,
к которому мы относим объект,
либо значение функции в задаче регрессии.                           

Объект -- набор вещественных чисел
(координаты обычно называют признаками, features).
В решающих деревьях рассматривают простейшие предикаты:
    Q(j, t):   x_j <= t
t <- R -- пороговое значение (threshold),
     j -- номер координаты объекта

S -- множество объектов, которые относятся к данной вершине
     дерева.
Предикат разбивает его на два подмножества     
     S_left  = { x <- S: Q(x_j, t) = 0, т.е. x_j > t }
     S_right = { x <- S: Q(x_j, t) = 1, т.е. x_j <= t }
     
Для построения решающего дерева используется жадный
алгоритм: у нас уже есть набор объектов, которые относятся
к данной вершине дерева. Мы хотим максимально хорошо разбить
наше множество на два подмножества так, чтобы разнообразие  
объектов внутри каждого подмножества максимально уменьшилось.
К сожалению, это приходится делать чистым перебором:
перебираются всевозможные пары (j, t), где j -- номер координаты
объектов, t -- пороговое значение для j-й координаты.
Для каждой пары вычисляем, насколько хорошее разбиение
(насколько сильно уменьшается разнообрази элементов
в каждом из подмножеств), и вбираем самое хорошее разбиение.     
     
                     Вершина дерева V
                     с множеством объектов S   
                     Выбрали предикат Q(j, t)
                       /         \
                   S_left =    S_right = { x<-S: Q(x_j, t) = 1 }
       = { x<-S: Q(x_j, t) = 0 }

Критерии остановки:
   -- все или почти все элементы
      принадлежат одному классу
      (или дисперсия значений в
      задаче регрессии невелика);
   -- |S| невелика
   -- высота дерева может быть ограничена
   
В задаче классификации для листа мы выбираем самый
популярный класс; в задаче регрессии выбирается
среднее значение функции для всех элементов,
входящих в лист.

Качество вершины дерева определяется разнообразием
входящих в нее элементов -- по-английски
     impurity (нечистота, грязь, наличие примесей)
     или по-русски информационное разнообразие,
                  информационный критерий.
                  
1. Рассмотрим сначала задачу регрессии.
     H(v) = 1/|S| Sum_(x,y)<-S L(y, c), где с -- значение ф-ции
                                        в вершине дерева,
                           L -- функция потери (Loss function)
Самая простая ф-ция потери -- квадрат разности
                 L(y, c) = (y - c)^2 
                 
Легко доказать, что функция H(v) минимальна при с, равном
среднему значению y для всех (x, y) <- S, тогда H(v)
есть среднее квадратичное отлонение от среднего значения,
т.е. дисперсия значений y для элементов (x, y) <- S.

2. Теперь рассмотрим задачу классификации.
Тогда вместо одного числа, которое мы приписываем
вершине дерева в задачи регрессии, мы рассматриваем
вертор вероятностей того, что элемент принадлежит
данному классу.
Пусть число классов равно K. Тогда значение вершину
     с = (c1, c2, ..., c_K), где c_i -- вероятность того,
          что объект, входящий в эту вершину, принадлежит
          классу i.  c1 + c2 + ... + c_K = 1
          
     c* = argmin_c 1/|S| Sum_(x,y) Sum_k (c_k - [y = k])^2           
                                                                         
Это значение минимально при
     с = (p1, p2, ..., p_K),   p_k = |{(x,y)<-S: y = k}|/|S|
                               доля элементов, принадлежащих 
                               классу k.
Критерий Gini (Джини):
     H(v) = Sum_k p_k*(1 - p_k)
     
Пусть c = (p1, p2, ..., p_K). Подставим эти значения в формулу
для impurity:                      
    1/|S| Sum_(x,y) Sum_k (p_k - [y = k])^2 =
    = 1/|S| Sum_k Sum_(x,y) (p_k - [y = k])^2 = 
Обозначим n_k -- число объектов, принадлежащих классу k
    = 1/|S| Sum_k (
                Sum_(x,y): y = k  (p_k - 1)^2  +               
                Sum_(x,y): y != k  p_k^2
            )
    = 1/|S| Sum_k (
                n_k*(p_k - 1)^2 + 
                (|S| - n_k) p_k^2
            )
     = Sum_k (
                (n_k*(p_k - 1)^2) / |S| +
                (|S| - n_k) p_k^2 / |S|
             )
     = Sum_k (
                p_k*(p_k - 1)^2 + 
                (1 - p_k)*p_k^2
             )
     = Sum_k (
                p_k^3 - 2*p_k^2 + p_k +
                p_k^2 - p_k^3
             )
     = Sum_k ( p_k - p_k^2) = Sum_k p_k*(1 - p_k)                                              

Получили критери Gini:
     H(v) = Sum_k p_k*(1 - p_k)
Он показывает, насколько велико разннобразие элементов
из нашей тренировочной выборки в вершине дерева (impurity).

Иногда вместо критерия Джини используется критерий энтропии.
Информация == энтропия:
есть распределение вероятностней p1, p2, ..., p_k
    H(p) = Sum p_i * (-log_2 p_i)
    
Откуда получается эта формула:
Пусть мы задаем вопрос и ответы да и нет равновероятны.
Тогда, получив ответ на вопрос, мы имеем 1 бит информации.
Пример: у нас есть 10 цифр 0, 1, ..., 9
Мы задаем два вопроса:
           число четное?
          да /      \ нет
          1/2     число меньше 5?
      0,2,4,6,  да  /   \ нет
      8           1/4   1/4 (чуть меньше)
              5,7,9      1,3

Количество информации 5*1 + 3*2 + 2*2 =
                      5*(-log_2(1/2)) + 
                      3*(-log_2(1/4)) +
                      2*(-log_2(1/4))
                      
Нормируем это выражение, получаем формулу для количества
информации (энтропии):
                      5/10*(-log_2(1/2)) + 
                      3/10*(-log_2(1/4)) +
                      2/10*(-log_2(1/4))
                = p1*(-log_2(p1)) + p2*(-log_2(p2)) +
                  + p3*(-log_2(p3))
                      
Критерий энтропии, можно использовать вместо критерия Джини:
    p_k -- доля объектов класса k, которые относятся
           к вершине дерева (всего S таких объектов)
    H(v) = -Sum_k p_k*log(p_k)           
    
Мы выбираем предикат Q(j, t) такой, что выражение
        H(v) - |left(v)| /|v| * H(left(v)) - 
             - |right(v)|/|v| * H(right(v))
максмально по всем парам (j, t).

Реализация на Питоне

from sklearn.tree import DecisionTreeClassifier
from sklearn.tree import export_text

        clf = DecisionTreeClassifier(max_depth = maxDepth)
        X = np.array([[p[0], p[1]] for p in points])
        y = np.array([c for c in targetClassNumbers])

        clf.fit(X, y)
        print(export_text(clf))

        y_predicted = clf.predict(X)
        
Раскрашивание плоскости:

      def onPaint():
        if clf == None:
            return
        dx = 0.2
        dy = 0.2
        x0 = xMin(); x1 = xMax()
        y = yMin()
        y1 = yMax()
        while y < y1:
            x = x0
            while x < x1:
                c = clf.predict([[x, y]])[0]
                color = alphaColor(c)
                leftTop = R2Point(x, y + dy)
                rightBottom = R2Point(x + dx, y)
                lt = map(leftTop)
                rb = map(rightBottom)
                rectID = drawArea.create_rectangle(
                    lt[0], lt[1], rb[0], rb[1],
                    fill = color, outline = color
                )
                rectIDs.append(rectID)
                x += dx
            y += dy
        drawPoints()
      
===================================================

12 May 2023

Решающие деревья (продолжение):
использование в задаче регрессии

Дата тренировочная база объектов
Объект -- вектор (набор признаков) n-мерного пространства.
Вектор обозначается буквой x.
В задаче регрессии мы строим функцию
    y = y(x)
    y: R^n --> R
В базе данных мы имеем набор пар
    (x_i, y_i) -- объект, значение функции.
Надо построить математическую модель, которая для произвольного
вектора x вычисляет функцию y = y(x).

Как выбирается предикат в вершине?
Наша цель -- уменьшить разнообразие элементов
(impurity) в двух дочерних вершинах.

Как измерить разнообразие?
В задаче регрессии это просто дисперсия значений y_i
для элементов (x_i, y_i), входящих в вершину v дерева.

В задаче классификации обычно применяется один из двух методов:
1) энтропия (количество информации)
   для каждого класса k вычисляется доля p_k элементов
   из класса k, принадлежащих этой вершине
   H(v) = Sum_k p_k*(-log_2 p_k);
2) критерий Джинни:
   H(v) = Sum_k p_k*(1 - p_k) = Sum_k (p_k - p_k^2) =
        = 1 - Sum_k (p_k)^2
        
В лототроне шарики классов 1, 2, ..., K,
номер класса внутри шарика (мы его не видим).
Мы выбираем случайный шарик и относим его к одну из классов
1, 2, ..., K в соответствии с вероятностями p_k, k = 1..K.
Какова вероятность ошибки? Как раз эта сумма:
    Sum_k p_k*(1 - p_k)
Или, иначе, вероятность правильного выбора равна
    Sum_k (p_k)^2  
Тогда вероятность ошибки
    1 - Sum_k (p_k)^2 = Sum_k p_k*(1 - p_k).       

В обеих задачах (регрессия и классификаци)
мы выбираем предикат в вершине v так, чтобы при переходе к
переходе к дочерним вершинам разнообразие максимально
уменьшилось:
                      v        
        Q(v)(x) = 1 /   \ Q(v)(x) = 0 
              left(v)  right(v)
Предикат Q(v) = Q(j_v, t_v)
         Q(v)(x) = (coord_j(x) <= t)

===========================================

19 May 2023

Схема градиентного бустинга в применении к решающим деревьям

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

Как можно бороться с переобучением деревьев?

Первая идея -- использовать несколько деревьев (лес, forest),
каждое из которых обучается на некотором подмножестве
нашей выборки. На самом деле, обычно применяют
    bootstrap -- выборка с возвращением.
    Пусть имеем L объектов. Выбираем слуайно первый объект,
    возвращаем его в базу, затем второй объект и т.д.
    Всего выбираем L объектов. При этом в построенной выборке
    некоторые объекты могут совпадать (например, мы выберем
    2/3 объектов из исходной выборки, некоторые объекты могут
    совпадать).
Всего создаем N подобных выборок X1, X2, ..., X_N,
|X_i| = L. По каждой выборке X_i строим модель b_i,
основаная модель будет суммой (вернее, средним значением)
построенных моделей:
    a(x) = 1/N (b1(x) + b2(x) + ... + b_N(x))
Смещение модели a будет таким же, как и у моделей b_i,
зато разброс моделей a в зависимости от случайных выборок будет
меньше: в идеальном случае, когда элементы выборок X_i, X_j
будут независимы друг от друга, разброс (дисперсия, переобучение)
уменьшится в N раз. Реальные случаи не идеальны, но для улучшения
независимости можно вычислять энтропию в родительской и дочерних
вершинах не на всём множестве объектов, а только на некотором 
подмножестве.

Однако, как оказалось, использование леса решающих деревьев --
не самый лучший способ. Существует другой метод, гораздо
более сильный: градиентный бустинг (Фридман, 1999).

1. Чтобы избежать переобучения, будем использовать
   простые деревья: либо
   1) максимальная глубина дерева невелика -- скажем,
      глубины 3 (по умолчанию), либо
   2) количество листьев не превышает заданного числа
      (скажем, 8),
   3) минимальное число объектов в листе (число объектов - 1,
      при котором мы еще можем разделить вершину дерева в
      две дочерние вершины), ...
      
2. Наша модель будет представлять собой ансамбль деревьев:
   a_N(x) = 1/N (b1(x) + b2(x) + ... + b_N(x))
При этом, если построена модель A_N-1(x), то мы к ней добавляем
простую модель b_N(x), которая строится так, чтобы она
исправила ошибки предудыщей (составной) модели A_N-1(x):
   Тренировочная выборка: (x_i, y_i), i = 1, 2, ..., L
   a_N(x_i) = y_i^
   Мы строим нашу модель b_N+1 по тренировочной базе
            (x_i, y_i - y_i^)
В случае функции ошибки MSE модель b_N+1 минимизирует
квадратичную функцию ошибки
   MSE(b) = 1/L Sum_i (b_N+1(x_i) - (y_i - y_i^))^2
В данном случае функция потери -- разность квадратов:
   L(y, z) = (y - z)^2
В общем случае функция потери (на одном элементе)
может быть другого вида, например, несимметричной.
   Error(b) = 1/L Sum_i L(b_N+1(x_i), (y_i - y_i^))        
                         
В случае деревьев алгоритм построения такого ансамбля деревьев
называется GBDT -- Gradient Boosting Decision Trees.