Практическое применение принципа максимума при математическом моделировании процессов. Принцип максимума л.С.Понтрягина в теории оптимальных систем

Учтем теперь ограничения (2.2.2) на управление. Если в процессе оптимального управления функции не достигают границ множества (2.2.2) (что означает ) то для них выполняются соотношения (2.2.13), (2.2.14). Однако часто оптимальное управление принимает граничные значения либо - , более того, оптимальное управление может скачком переходить с одной границы на другую. Такие управления уже являются кусочно-непрерывными функциями времени.

При попадании оптимального управления на границу множества U соотношения (2.2.13), (2.2.14) нарушаются. Оптимальные управления удовлетворяют в этом случае принципу максимума Л. С. Понтрягина, установленного и доказанного в форме приведенной ниже - теоремы.

Переходя к этой теореме, сделаем некоторые пояснения. Возьмем произвольное допустимое управление и при начальных условиях найдем решение системы (2.2.1): .

Подставляя это решение и управление в (2.2.8), определим, пока при некоторых произвольных начальных условиях , решение (2.2.8): . При фиксированных (постоянных) значениях векторов функция Н становится функцией вектора . Максимум этой функции по и обозначим через :

Максимум (наибольшее значение) непрерывной функции может достигаться как в точках локального максимума этой функции, в которых

так и на границах множества .

Теорема 2.2.1 (принцип максимума Л. С. Понтрягина). Пусть , - такое допустимое управление, что соответствующие ему решения уравнения (2.2.11), исходящие в момент из состояния (2.2.3), (2.2.7), проходят в момент времени через точку .

Для оптимальности управления (при котором ) принимает наименьшее значение) необходимо существование таких ненулевых непрерывных функций удовлетворяющих уравнениям (2.2.12), что при любом функция переменного достигает при максимума

при этом в конечный момент времени выполняются соотношения

Если удовлетворяют (2.2.11), (2.2.12) и (2.2.17), то функции переменного t являются постоянными и поэтому проверку соотношений (2.2.18) можно проводить не обязательно в момент времени а в любой момент .

Доказательство теоремы является достаточно сложным, и поэтому в приложении 2 приведен лишь вывод основного соотношения (2.2.17) теоремы для случая свободного правого конца ( не задан) и фиксированного .

Соотношения (2.2.17) и (2.2.18) можно записать в более простой форме:

Таким образом, центральным в теореме 2.2.1 является условие максимума (2.2.19). Оно означает, что если - оптимальные управления, а -оптимальные траектории, то непременно найдутся такая постоянная и такие решения ), системы (2.2.12), что функция их, переменных при всех будет достигать максимума на U именно при оптимальных управлениях . Поэтому теорему 2.2.1, дающую необходимое условие оптимальности в задачах оптимального управления, принято называть принципом максимума. Отметим, что во внутренних точках множества U для оптимального управления выполняются условия (2.2.13), (2.2.14), которые являются необходимыми для (2.2.19).

Практическое применение принципа максимума.

Как же практически воспользоваться условием (2.2.19), ведь функции и постоянная , входящие в это условие, неизвестны?

Здесь поступают следующим образом: рассматривая функцию ) как функцию переменных и считая переменные параметрами, решают задачу максимизации функции и находят функцию

на которой достигается наибольшее значение функции .

В ряде случаев функция (2.2.20) может быть записана в явном виде. Например, если правые части (2.2.1) имеют структуру

а подынтегральное выражение функционала (2.2.5)

множество описывается U неравенствами (2.2.2), то

и эта функция достигает наибольшего значения на U в точке с координатами

Формула (2.2.22) дает большой объем информации о структуре оптимального управления: координата оптимального управления является ступенчатой (кусочно-постоянной) функцией со значениями при этом моменты переключения определяются условием

Итак, допустим, что функция (2.2.20) известна. Рассмотрим систему дифференциальных уравнений

Функции и , входящие в правые части этих уравнений, известны. Общее решение системы (2.2.24), (2.2.25) зависит от произвольных постоянных, которые определяются из краевых условий (2.2.3), (2.2.4). Задача интегрирования уравнений (2.2.24), (2.2.25) при краевых условиях (2.2.3), (2.2.4) называется краевой задачей (двухточечной краевой задачей).

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

Трудность ее решения состоит в том, что интегрирование уравнений (2.2.24), (2.2.25) в «прямом времени» не представляется возможным, так как неизвестны начальные условия Один из возможных подходов к решению краевой задачи заключается в следующем. Задаваясь произвольным вектором и интегрируя (2.2.24), (2.2.25) при известных начальных условиях , найдем функции и при проверим выполнение равенства (2.2.4). Если оно нарушается, задаемся другим вектором и, интегрируя (2.2.24), (2.2.25) при начальных условиях , получим при вектор .

Если он не совпадает с заданным, продолжаем процесс до тех пор, пока не найдется такой вектор что условия (2.2.4) будут выполняться с приемлемой точностью. При этом подходе используются градиентные методы, когда определяется из условия минимума «расстояния» от заданного вектора .

В вычислительной математике разработан ряд методов приближенного численного решения краевых задач: метод стрельбы, метод прогонки, ряд итерационных методов , . Во многих случаях не представляется возможным найти из условия (2.2.19) явный вид (2.2.22) оптимального управления. Тогда уравнения (2.2.1), (2.2.6), сопряженная система (2.2.12) и условия максимума (2.2.19) образуют краевую задачу принципа максимума. Эта задача имеет ряд специфических особенностей, затрудняющих применение стандартных численных методов решения краевых задач. К числу таких особенностей относятся разрывы функций удовлетворяющих условию максимума (2.2.14), их неединственность, нелинейный характер зависимости (2.2.20) даже в линейных системах. Кроме того, особенностью краевых задач, связанных с принципом максимума даже в случаях, когда удается найти явный вид управлений (2.2.20), является их плохая сходимость, вызванная неустойчивостью системы (2.2.24), (2.2.25). Ряд приемов решения краевых задач принципа максимума изложен, например, в .

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

Пример 2.2.1. Построение оптимального по расходу топлива управления .

Рассмотрим объект управления, описываемый уравнениями

Пусть на управление наложено ограничение

Функционал оптимизации, выражающий расход топлива, имеет вид

Заданы начальное состояние

и условие в момент времени

Требуется найти , при котором объект (2.2.26) переходит из состояния (2.2.29) в состояние (2.2.30), при этом выполняются ограничения (2.2.27), а функционал (2.2.28) принимает наименьшее значение.

Переходя к определению оптимального управления на основе принципа максимума, сформируем функцию

уравнения для вспомогательных переменных

Управление , доставляющее максимум функции (2.2.31), определяется как

Уравнения (2.2.26), (2.2.32), (2.2.33) составляют краевую задачу. Переходя к ее исследованию, запишем решение системы (2.2.32):

где - неизвестные числа, которые необходимо определить так, чтобы управление (2.2.33) привело объект (2.2.26) в состояние (2.2.30).

Найдем решение системы (2.2.26) при и . В первом случае решение этой системы имеет вид . Оно зависит от постоянных R и р, при этом . Фазовые траектории этой системы представляют собой окружности с центром в начале координат (рис. 2.2.1, а). Фазовые траектории системы (2.2.26) при и также являются окружностями, центры которых расположены в точках соответственно (рис. 2.2.1, б, в).

В этом параграфе мы докажем свойство решений одномерного уравнения теплопроводности, которое называется принципом максимального значения . Оно может быть сформулировано как теорема.

Т е о р е м а. Если функция u (x ,t ), определенная и непрерывная в замкнутой области и , удовлетворяет в этой области уравнению теплопроводности

то максимальное и минимальное значения функции u (x ,t ) достигаются или в начальный момент времени или в граничных точках x = 0 или x = l.

Функция , очевидно, удовлетворяет уравнению (40) и достигает своего максимального (минимального) значения в любой точке. Однако это не противоречит теореме, так как из её условия следует, что если максимальное (минимальное) значение достигается внутри области, то оно также должно достигаться или при t= 0, или при x = 0 илипри x=l.

Физический смысл этой теоремы очевиден и заключается в следующем. Если температура на границе или в начальный момент не превосходит некоторого значения M , то при отсутствии источников тепла внутри тела не может создаться температура, больше чем М .

Остановимся на доказательстве теоремы для максимального значения. Оно ведется от противного. Итак, пусть М максимальное значение функции u (x ,t ) при t = 0 (0 ≤ x l ) или при x = 0 илипри x = l (0 ≤ t T ). Допустим теперь, что в некоторой точке области (x 0 , t 0), такой, что 0 < x 0 < l и0 < t 0 ≤ T , функция u (x ,t ) достигает своего максимального значения, превосходящего М на величину ε, т.е.

Тогда в точке (x 0 , t 0) должны выполняться соотношения

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

где k – постоянный коэффициент. Очевидно, что

Выберем так, чтобы kT было меньше ε/2, т.е. , тогда максимальное значение v (x , t ) при t = 0 (0 ≤ x l ) или при x = 0 илипри x = l не будет превосходить , т.е.

(при t = 0 или x = 0 или x = l ), (44)

так как для этих аргументов первое слагаемое в формуле (43) не превосходит М , а второе .

В силу непрерывности функции v (x , t ), она должна в некоторой точке (x 1 ,t 1) достигать своего максимального значения, причем

Момент времени t 1 строго больше нуля и , так как при или , или имеет место неравенство (44). В точке (x 1 , t 1), по аналогии с (41) и (42), должно быть

Имея в виду определение функции v (x , t ) (43), получим

Отсюда следует, что

т.е. уравнение (40) во внутренней точке (x 1 ,t 1) не удовлетворяется. Тем самым доказано, что решение u (x ,t ) уравнения теплопроводности (40) внутри области не может принимать значений, превосходящих наибольшее значение u (x ,t ) на границе.



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

Приведем и докажем следствия из принципа максимального значения:

Следствие 1. Если два решения уравнения (40) и удовлетворяют условиям:

,

Доказательство. В силу линейности (40) функция является его решением, следовательно, удовлетворяет принципу максимального значения. При этом:

Следовательно:

в противном случае имела бы отрицательное минимальное значение. Следствие 1 доказано.

Следствие 2. Если три решения уравнения (40) , и удовлетворяют условию:

при , и , то это же неравенство выполняются и для всех .

Доказательство. Проводится просто применением следствия 1 к парам функций и , и .

Рассмотрим в решение уравнения (1), соответствующее начальному и граничным условиям вида:

Пусть есть решение уравнения (40), соответствующее возмущенным начальному и граничным условиям, задаваемыми функциями , и , такими, что:

Используем следствие 3, можем заключить, что: , что и подразумевает сколь угодную близость решений исходной и возмущенной задач.

Пример

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

  • X количество успехов в двенадцати независимых испытаний с распределением Бернулли с вероятностью успеха θ в каждом из них.
  • Y количество независимых испытаний с распределением Бернулли, необходимых для получения трех успехов. Вероятность успеха в каждом из испытаний θ.

Тогда рассмотрение X = 3 даст функцию правдоподобия

а рассмотрение Y = 12 даст функцию правдоподобия

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

Разница в наблюдении X = 3 и наблюдении Y = 12 исключительно в дизайне эксперимента: в одном случае изначально было решено делать двенадцать попыток, а в другом делать попытки, пока не будет трех успешных. Результат будет одинаковым в обоих случаях. Поэтому принцип максимального правдоподобия иногда выражают следующим образом:

Вывод должен зависеть только от исхода эксперимента, а не от дизайна эксперимента.

Закон максимального правдоподобия

Связанная с принципом максимального правдоподобия концепция - это закон максимального правдоподобия , говорящий, что отношение того, какое значение параметра более применимо, равняется отношению их функций правдоподобия. Тогда отношение

является мерой того, насколько величина x принимает параметр a в отношении к b . Таким образом, если отношение равняется 1, то разницы нет, а если больше 1, то a предпочтительней b , и наоборот.

Из принципа максимального правдоподобия и закона максимального правдоподобия следует, что параметр, который максимизирует функцию правдоподобия, является лучшим. Это и является основой широко известного метода максимального правдоподобия .

Историческая справка

Принцип максимального правдоподобия был впервые упомянут в печати в г. Однако основы принципа и применение его на практике были опубликованы ранее в работах Р. А. Фишера в г.

Аргументы за и против принципа максимального правдоподобия

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

Зависимость результата от организации эксперимента

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

Некоторые классичекие методы проверки гипотез базируются не на правдоподобии. Часто приводимый пример это проблема оптимальной остановки. Предположим я сказал, что бросил монету 12 раз и получил 3 решки. Из этого вы сможете сделать некоторые выводы о вероятности выпадения решки у этой монеты. А теперь предположим, что я бросал монету пока решка не выпала 3 раза, и получилось 12 бросков. Сделаете ли вы теперь другие выводы?

Функция правдоподобия одинакова в обоих случаях и пропорциональна

.

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

Предположим некоторая группа ученых определяет вероятность некоторого исхода (который мы будем называть "успехом") серией экспериментов. Здравый смысл подсказывает нам, что если нет оснований считать что успех более вероятен, чем неудача, и наоборот, то следует положить вероятность успеха равной 0.5. Ученый Адам сделал 12 испытаний, в которых получил 3 успеха и 9 неудач, после чего умер.

Его коллега по лаборатории Билл продолжил работу Адама и опубликовал результат проверки гипотезы. Он проверил гипотезу что вероятность успеха p =0.5 против p < 0.5. Вероятность того, что в 12 испытаниях наступит не более 3 успехов, равна

что есть 299/4096 = 7.3 %. Таким образом гипотеза не отвергается при 5 % уровне доверия.

Шарлотта, прочитав статью Билла, пишет письмо. Она считает, что Адам, возможно, продолжал испытания пока не умер, успев получить к этому моменту 3 успеха. Вероятность того, что для трех успехов потребуется 12 или более испытаний равна

что есть 134/4096 = 3.27 %. И теперь результат отвергается при уровне в 5 %.

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

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

Литература

См. также

Ссылки

  • Anthony W.F. Edwards. «Likelihood». http://www.cimat.mx/reportes/enlinea/D-99-10.html
  • Jeff Miller. Earliest Known Uses of Some of the Words of Mathematics (L)
  • John Aldrich. Likelihood and Probability in R. A. Fisher’s Statistical Methods for Research Workers

Wikimedia Foundation . 2010 .

Постановка задачи. Пусть состояние управляемой системы характеризуется п-мерным вектором состояния x(t). Целенаправленное воздействие на процесс можно осуществлять с помощью m-мерного вектора управления u(t). На векторы управления и состояния могут накладываться ограничения: x(t)eX, u(t)eU, где U и х - соответственно области допустимых управлений и состояний. Будем считать допустимыми управляющими воздействиями кусочно-непрерывные функции на отрезке управления , в точках разрыва первого рода непрерывные справа и на концах отрезка. Между x(t) и u(t) существует зависимость, записываемая в виде системы дифференциальных уравнений

Заданы граничные условия: начальное состояние управляемой системы x,(t 0) = *f и конечное состояние -х { (??) = **, i = l,2,...,n в n-мерном фазовом пространстве. Момент времени п полагаем незафиксированным, а характеризующим только момент перехода в конечное состояние х 1 . Критерий качества течения процесса управления описывается функционалом

Теперь постановку задачи оптимального управления можно сформулировать следующим образом. В п-мерном фазовом пространстве X даны две точки, х° и х 1 . Среди всех допустимых управлений, для которых фазовая траектория, исходящая в момент времени t 0 из точки х°, приходит в точку х 1 , найти такое управление, для которого функционал I принимает экстремальное значение. Такое управление называется оптимальным управлением, а соответствующая ему траектория - оптимальной траекторией.

Постановку задачи можно видоизменить, введя в рассмотрение еще одну координату вектора состояния, характеризующую текущее значение функционала:

Дифференцируя, получаем уравнение относительно новой координаты вектора состояния

I

с граничными условиями

Теперь речь будет идти о (п + 1)-мерном фазовом пространстве состояния X, а постановку задачи, эквивалентную предыдущей, можно сформулировать следующим образом.

В (п + 1)-мерном фазовом пространстве X заданы: точка с координатами (0,х°) и прямая П, проходящая параллельно оси х 0 через точку (0,х х). Среди всех допустимых управлений, для которых соответствующая фазовая траектория, исходя в момент времени из точки (0,х°), пересечет прямую П, найти управление, обеспечивающее наименьшее возможное значение координаты пересечения с прямой П вдоль оси х 0 .

Постановка задачи геометрически интерпретирована для случая п = 2 на рис. 2.3.

Кроме основной системы дифференциальных уравнений (2.22) и (2.24), которые вместе запишутся как

введем в рассмотрение дополнительную систему дифференциальных уравнений относительно вспомогательной вектор-функции |/(t)


Рис. 2.3.

После введения понятия функции Гамильтона

системы (2.25) и (2.26) можно объединить записью в виде так называемой гамильтоновой системы

Обозначим через М(х, i) максимальное значение функции Гамильтона

или, точнее, значение ее верхней грани - супремума.

Основная теорема. Пусть u(t) на отрезке - допустимое управление, удовлетворяющее постановке задачи. Тогда для оптимальности вектора управления u(t) необходимо, чтобы существовала ненулевая вектор-функция |/(t), такая, что:

1) для всех t на отрезке ] функция Гамильтона как функция и, ueU достигала максимума

2) в конечный момент времени t = П выполнялись соотношения

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

Таким образом, мы имеем 2(п + 1) + ш соотношений (2.25), (2.26) и (2.29), между 2(п + 1) + ш координатами векторов x(t), p(t) и u(t). Так как т соотношений не дифференциальные, то решение системы (2.25), (2.26) и (2.29) зависит от 2(n + 1) неизвестных параметров, кроме того, >t 0 или ti -t 0 = т также является параметром. Один из параметров несуществен, так как v|/(t) определяется с точностью до общего множителя в силу однородности Я относительно |/. Таким образом, для нахождения 2(n + 1) параметров имеем (2n + 1) граничных условий на функцию x(t) и второй пункт принципа максимума.

Принцип максимума для оптимального быстродействия. Важным для технических приложений является класс задач на оптимальное быстродействие: требуется найти управление u(t), переводящее фазовую точку из состояния х° в момент времени t 0 в х 1 за минимальное время, т.е. функционал в данном случае запишется как

т.е. / 0 (x,u)sl. Функция Гамильтона принимает вид

где

Таким образом, величина р 0 не влияет на значение u(t), при котором достигается максимум, а влияет только на величину максимума. В силу сказанного в принципе максимума для оптимального быстродействия можно говорить о функциях Н(х, ц/, и) и М(х, |/, и) = max Н(х, |/, и), поэтому равенство (2.29) основной теоремы будет иметь вид

а равенство (2.30) будет иметь вид

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

Полагаем, что управляющие воздействия подчинены ограничениям вида

Требуется минимизировать время перехода из х° в х 1 , т.е. I = t x -t 0 . Составим функцию Гамильтона:

Изменим порядок суммирования во второй группе слагаемых:

На основании принципа максимума u(t) следует искать, исходя из максимума Н(х, |/, и) относительно u (О для всех f из при учете ограничения (2.31). Нетрудно видеть (рис. 2.4), что значение, прини-

маемое u ; (t) в каждый момент t, определяется знаком?b^|/,-(t), т.е. изменение u, (t) происходит по закону,=1

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

Изменим порядок суммирования в первой группе слагаемых выражения (2.32)

Уравнения гамильтоновой системы относительно вспомогательной функции j/(t) (2.26) будут иметь вид


Рис. 2.4.

t*, t** - точки переключения управления

Система уравнений является однородной, следовательно, общее решение можно записать в виде

где [Ху] - совокупность корней характеристического уравнения или собственные значения матрицы А = [а,-у ].

Полагаем, что Xj,j=l,2,...,n являются простыми вещественными корнями. Тогда каждая из функций v|/, (t) как сумма п монотонных функций не более (п - 1) раза пересекает ось t. Так как функция

Cy(t) = ^byyl|/y(t) является суммой п монотонных функций, то можно 1=1

сформулировать второе свойство задачи на линейное оптимальное быстродействие систем n-го порядка, корни характеристического уравнения которых вещественны: управляющие воздействия имеют не более п промежутков знакопостоянства или не более (п - 1) переключения.

Поясним результат на конкретном примере. Пусть управляемая система представляет собой двойное интегрирующее звено, т.е.

с ограничением на управляющее воздействие | u(t) |

Требуется найти u(t), переводящее фазовую точку из произвольного положения х° в начало координат фазового пространства за минимальное время.

1. Введем фазовые координаты:

Тогда в нормальной форме уравнения системы запишутся как

2. Запишем выражение для функции Гамильтона:

Гамильтонова система относительно вспомогательной вектор-функции j/(t) имеет вид

Отсюда

3. На основании свойств управления:

Уравнения семейств фазовых траекторий (рис. 2.5, а, б):

Под действием управления н = ±1 фазовая точка может попасть в начало координат только по выделенной траектории (рис. 2.6). Чтобы фазовая точка попала в начало координат не более чем за одно переключение, движение должно быть организовано, как показано на рис. 2.6.

Задача оптимального управления конечным состоянием. Задача управления конечным состоянием при отсутствии ограничений на вектор управления - это задача Майера классического вариационного исчисления. В этих задачах функционал определяет значение функции от координат состояния в момент времени t = t b т.е.


Рис. 2.5.

а - u(t) = + 1; б - u(t) = -1


Рис. 2.6.

Пусть управляемая система описывается системой дифференциальных уравнений n-го порядка

с граничными условиями

Эти уравнения линейны относительно m-мерного вектора управления u(t). На управляющие воздействия накладываются ограничения вида

Требуется определить вектор управления, минимизирующий функционал

Введем новую координату

Дополнительное дифференциальное уравнение относительно x 0 (t) запишется в виде

с граничными условиями

Функция Гамильтона для данной системы

Выполним следующие преобразования:

При выполнении преобразований в силу ф; = 0 и |/ 0 (П) = -1 полагаем v|/ 0 (t) = -l. На основании принципа максимума при указанных ограничениях на управление оптимальное управление конечным состоянием для данной системы определяется выражением

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

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

Задача управления на минимум расхода энергии. Эта задача оптимального управления имеет большое практическое значение. Расход энергии, затраченной на управление, пропорционален интегралу по времени от квадрата управляющего воздействия. Если расход энергии по всем входам брать с одинаковыми весами, то функционал можно записать в виде

где коэффициент 1/2 введен для удобства последующих выкладок.

Рассмотрим систему, описываемую системой дифференциальных уравнений п-го порядка:

Рассмотрим два случая определения оптимального управления по расходу энергии на управление:

1) когда на вектор управления не накладывается ограничений;

2) когда управляющие воздействия подчиняются ограничениям вида

Введем новую координату

Соответствующее этой координате дифференциальное уравнение при граничных условиях

Функция Гамильтона для этой системы имеет вид

Меняем порядок суммирования в последней группе слагаемых и полагаем, как и прежде, j/ 0 =-1. Тогда

1-й случай. На вектор u(t) не накладывается ограничений.

Максимум функции Гамильтона Н{х, |/, и ) относительно и может быть найден на основании необходимого условия экстремума функции из классического анализа

Дифференцируем по п; :

Отсюда оптимальное управление

2-й случай. Управляющие воздействия подчиняются ограничению

Из предыдущего случая ясно, что, если | С; (?) | Uj(t) = С другой стороны, для тех t, для которых | Cj(t ) | > Mj, с учетом ограничения на Uj(t) функция Гамильтона максимальна, если Uj(t ) = MjSignCj(t). Отсюда заключаем, что оптимальное управление может быть записано в виде

Для определения Cj(t) необходимо решить гамильтонову систему уравнений.

Метод наименьших квадратов

Рассказ про линейные модели мы начнем с линейной регрессии. В первую очередь, необходимо задать модель зависимости объясняемой переменной от объясняющих ее факторов, функция зависимости будет линейной: . Если мы добавим фиктивную размерность для каждого наблюдения, тогда линейную форму можно переписать чуть более компактно, записав свободный член под сумму: . Если рассматривать матрицу наблюдения-признаки, у которой в строках находятся примеры из набора данных, то нам необходимо добавить единичную колонку слева. Зададим модель следующим образом:





Можем выписать выражение для каждого конкретного наблюдения



Также на модель накладываются следующие ограничения (иначе это будет какая то другая регрессия, но точно не линейная):



Оценка весов называется линейной, если



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



Один из способов вычислить значения параметров модели является метод наименьших квадратов (МНК), который минимизирует среднеквадратичную ошибку между реальным значением зависимой переменной и прогнозом, выданным моделью:



Для решения данной оптимизационной задачи необходимо вычислить производные по параметрам модели, приравнять их к нулю и решить полученные уравнения относительно (матричное дифференцирование неподготовленному читателю может показаться затруднительным, попробуйте расписать все через суммы, чтобы убедиться в ответе):


Шпаргалка по матричным производным




Итак, имея в виду все определения и условия описанные выше, мы можем утверждать, опираясь на теорему Маркова-Гаусса , что оценка МНК является лучшей оценкой параметров модели, среди всех линейных и несмещенных оценок, то есть обладающей наименьшей дисперсией.

Метод максимального правдоподобия

У читателя вполне резонно могли возникнуть вопросы: например, почему мы минимизируем среднеквадратичную ошибку, а не что-то другое. Ведь можно минимизировать среднее абсолютное значение невязки или еще что-то. Единственное, что произойдёт в случае изменения минимизируемого значения, так это то, что мы выйдем из условий теоремы Маркова-Гаусса, и наши оценки перестанут быть лучшими среди линейных и несмещенных.


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


Как-то после школы я заметил, что все помнят формулу этилового спирта. Тогда я решил провести эксперимент: помнят ли люди более простую формулу метилового спирта: . Мы опросили 400 человек и оказалось, что формулу помнят всего 117 человек. Разумно предположить, что вероятность того, что следующий опрошенный знает формулу метилового спирта – . Покажем, что такая интуитивно понятная оценка не просто хороша, а еще и является оценкой максимального правдоподобия.


Разберемся, откуда берется эта оценка, а для этого вспомним определение распределения Бернулли : случайная величина имеет распределение Бернулли, если она принимает всего два значения ( и с вероятностями и соответственно) и имеет следующую функцию распределения вероятности:



Похоже, это распределение – то, что нам нужно, а параметр распределения и есть та оценка вероятности того, что человек знает формулу метилового спирта. Мы проделали независимых экспериментов, обозначим их исходы как . Запишем правдоподобие наших данных (наблюдений), то есть вероятность наблюдать 117 реализаций случайной величины и 283 реализации :






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




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



но будем теперь считать, что случайные ошибки берутся из центрированного нормального распределения :



Перепишем модель в новом свете:



Так как примеры берутся независимо (ошибки не скоррелированы – одно из условий теоремы Маркова-Гаусса), то полное правдоподобие данных будет выглядеть как произведение функций плотности . Рассмотрим логарифм правдоподобия, что позволит нам перейти от произведения к сумме:



Мы хотим найти гипотезу максимального правдоподобия, т.е. нам нужно максимизировать выражение , а это то же самое, что и максимизация его логарифма. Обратите внимание, что при максимизации функции по какому-то параметру можно выкинуть все члены, не зависящие от этого параметра:



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

Разложение ошибки на смещение и разброс (Bias-variance decomposition)

Поговорим немного о свойствах ошибки прогноза линейной регрессии (в принципе эти рассуждения верны для всех алгоритмов машинного обучения). В свете предыдущего пункта мы выяснили, что:



Тогда ошибка в точке раскладывается следующим образом:



Для наглядности опустим обозначение аргумента функций. Рассмотрим каждый член в отдельности, первые два расписываются легко по формуле :



Пояснения:




И теперь последний член суммы. Мы помним, что ошибка и целевая переменная независимы друг от друга:



Наконец, собираем все вместе:



Итак, мы достигли цели всех вычислений, описанных выше, последняя формула говорит нам, что ошибка прогноза любой модели вида складывается из:



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


Как правило, при увеличении сложности модели (например, при увеличении количества свободных параметров) увеличивается дисперсия (разброс) оценки, но уменьшается смещение. Из-за того что тренировочный набор данных полностью запоминается вместо обобщения, небольшие изменения приводят к неожиданным результатам (переобучение). Если же модель слабая, то она не в состоянии выучить закономерность, в результате выучивается что-то другое, смещенное относительно правильного решения.


Теорема Маркова-Гаусса как раз утверждает, что МНК-оценка параметров линейной модели является самой лучшей в классе несмещенных линейных оценок, то есть с наименьшей дисперсией. Это значит, что если существует какая-либо другая несмещенная модель тоже из класса линейных моделей, то мы можем быть уверены, что .

Регуляризация линейной регрессии

Иногда бывают ситуации, когда мы намеренно увеличиваем смещенность модели ради ее стабильности, т.е. ради уменьшения дисперсии модели . Одним из условий теоремы Маркова-Гаусса является полный столбцовый ранг матрицы . В противном случае решение МНК не существует, т.к. не будет существовать обратная матрица Другими словами, матрица будет сингулярна, или вырожденна. Такая задача называется некорректно поставленной . Задачу нужно скорректировать, а именно, сделать матрицу невырожденной, или регулярной (именно поэтому этот процесс называется регуляризацией). Чаще в данных мы можем наблюдать так называемую мультиколлинеарность - когда два или несколько признаков сильно коррелированы, в матрице это проявляется в виде "почти" линейной зависимости столбцов. Например, в задаче прогнозирования цены квартиры по ее параметрам "почти" линейная зависимость будет у признаков "площадь с учетом балкона" и "площадь без учета балкона". Формально для таких данных матрица будет обратима, но из-за мультиколлинеарности у матрицы некоторые собственные значения будут близки к нулю, а в обратной матрице появятся экстремально большие собственные значения, т.к. собственные значения обратной матрицы – это . Итогом такого шатания собственных значений станет нестабильная оценка параметров модели, т.е. добавление нового наблюдения в набор тренировочных данных приведёт к совершенно другому решению. Иллюстрации роста коэффициентов вы найдете в одном из наших прошлых постов . Одним из способов регуляризации является регуляризация Тихонова , которая в общем виде выглядит как добавление нового члена к среднеквадратичной ошибке:



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



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


Такое решение уменьшает дисперсию, но становится смещенным, т.к. минимизируется также и норма вектора параметров, что заставляет решение сдвигаться в сторону нуля. На рисунке ниже на пересечении белых пунктирных линий находится МНК-решение. Голубыми точками обозначены различные решения гребневой регрессии. Видно, что при увеличении параметра регуляризации решение сдвигается в сторону нуля.


Советуем обратиться в наш прошлый пост за примером того, как регуляризация справляется с проблемой мультиколлинеарности, а также чтобы освежить в памяти еще несколько интерпретаций регуляризации.

2. Логистическая регрессия

Линейный классификатор

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



Мы уже знакомы с линейной регрессией и методом наименьших квадратов. Рассмотрим задачу бинарной классификации, причем метки целевого класса обозначим "+1" (положительные примеры) и "-1" (отрицательные примеры).
Один из самых простых линейных классификаторов получается на основе регрессии вот таким образом:




Логистическая регрессия как линейный классификатор

Логистическая регрессия является частным случаем линейного классификатора, но она обладает хорошим "умением" – прогнозировать вероятность отнесения примера к классу "+":



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



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


Итак, мы хотим прогнозировать вероятность , а пока умеем строить линейный прогноз с помощью МНК: . Каким образом преобразовать полученное значение в вероятность, пределы которой – ? Очевидно, для этого нужна некоторая функция В модели логистической регрессии для этого берется конкретная функция: . И сейчас разберемся, каковы для этого предпосылки.



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


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


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




В правой части мы получили как раз сигмоид-функцию.


Итак, логистическая регрессия прогнозирует вероятность отнесения примера к классу "+" (при условии, что мы знаем его признаки и веса модели) как сигмоид-преобразование линейной комбинации вектора весов модели и вектора признаков примера:



Следующий вопрос: как модель обучается? Тут мы опять обращаемся к принципу максимального правдоподобия.

Принцип максимального правдоподобия и логистическая регрессия

Теперь посмотрим, как из принципа максимального правдоподобия получается оптимизационная задача, которую решает логистическая регрессия, а именно, – минимизация логистической функции потерь.
Только что мы увидели, что логистическая регрессия моделирует вероятность отнесения примера к классу "+" как



Тогда для класса "-" аналогичная вероятность:



Оба этих выражения можно ловко объединить в одно (следите за моими руками – не обманывают ли вас):



Выражение называется отступом (margin ) классификации на объекте (не путать с зазором (тоже margin), про который чаще всего говорят в контексте SVM). Если он неотрицателен, модель не ошибается на объекте , если же отрицателен – значит, класс для спрогнозирован неправильно.
Заметим, что отступ определен для объектов именно обучающей выборки, для которых известны реальные метки целевого класса .


Чтобы понять, почему это мы сделали такие выводы, обратимся к геометрической интерпретации линейного классификатора. Подробно про это можно почитать в материалах Евгения Соколова.



Ответ



Когда получим (или посмотрим) ответ, то поймем, что чем больше по модулю выражение , тем дальше точка находится от плоскости


Значит, выражение – это своего рода "уверенность" модели в классификации объекта :




Теперь распишем правдоподобие выборки, а именно, вероятность наблюдать данный вектор у выборки . Делаем сильное предположение: объекты приходят независимо, из одного распределения (i.i.d. ). Тогда



где – длина выборки (число строк).


Как водится, возьмем логарифм данного выражения (сумму оптимизировать намного проще, чем произведение):



То есть в даном случае принцип максимизации правдоподобия приводит к минимизации выражения



Это логистическая функция потерь, просуммированная по всем объектам обучающей выборки.


Посмотрим на новую фунцию как на функцию от отступа: . Нарисуем ее график, а также график 1/0 функциий потерь (zero-one loss ), которая просто штрафует модель на 1 за ошибку на каждом объекте (отступ отрицательный): .



Картинка отражает общую идею, что в задаче классификации, не умея напрямую минимизировать число ошибок (по крайней мере, градиентными методами это не сделать – производная 1/0 функциий потерь в нуле обращается в бесконечность), мы минимизируем некоторую ее верхнюю оценку. В данном случае это логистическая функция потерь (где логарифм двоичный, но это не принципиально), и справедливо



где – попросту число ошибок логистической регрессии с весами на выборке .


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

-регуляризация логистических потерь

L2-регуляризация логистической регрессии устроена почти так же, как и в случае с гребневой (Ridge регрессией). Вместо функционала минимизируется следующий:



В случае логистической регрессии принято введение обратного коэффициента регуляризации . И тогда решением задачи будет



3. Наглядный пример регуляризации логистической регрессии

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


Посмотрим, как регуляризация влияет на качество классификации на наборе данных по тестированию микрочипов из курса Andrew Ng по машинному обучению.
Будем использовать логистическую регрессию с полиномиальными признаками и варьировать параметр регуляризации C.
Сначала посмотрим, как регуляризация влияет на разделяющую границу классификатора, интуитивно распознаем переобучение и недообучение.
Потом численно установим близкий к оптимальному параметр регуляризации с помощью кросс-валидации (cross-validation) и перебора по сетке (GridSearch).


Подключение библиотек

from __future__ import division, print_function # отключим всякие предупреждения Anaconda import warnings warnings.filterwarnings("ignore") %matplotlib inline from matplotlib import pyplot as plt import seaborn as sns import numpy as np import pandas as pd from sklearn.preprocessing import PolynomialFeatures from sklearn.linear_model import LogisticRegression, LogisticRegressionCV from sklearn.model_selection import cross_val_score, StratifiedKFold from sklearn.model_selection import GridSearchCV


Загружаем данные с помощью метода read_csv библиотеки pandas. В этом наборе данных для 118 микрочипов (объекты) указаны результаты двух тестов по контролю качества (два числовых признака) и сказано, пустили ли микрочип в производство. Признаки уже центрированы, то есть из всех значений вычтены средние по столбцам. Таким образом, "среднему" микрочипу соответствуют нулевые значения результатов тестов.

data = pd.read_csv("../../data/microchip_tests.txt", header=None, names = ("test1","test2","released")) # информация о наборе данных data.info()


RangeIndex: 118 entries, 0 to 117
Data columns (total 3 columns):
test1 118 non-null float64
test2 118 non-null float64
released 118 non-null int64
dtypes: float64(2), int64(1)
memory usage: 2.8 KB


Посмотрим на первые и последние 5 строк.



Сохраним обучающую выборку и метки целевого класса в отдельных массивах NumPy. Отобразим данные. Красный цвет соответствует бракованным чипам, зеленый – нормальным.


Код

X = data.ix[:,:2].values y = data.ix[:,2].values
plt.scatter(X, X, c="green", label="Выпущен") plt.scatter(X, X, c="red", label="Бракован") plt.xlabel("Тест 1") plt.ylabel("Тест 2") plt.title("2 теста микрочипов") plt.legend();


Определяем функцию для отображения разделяющей кривой классификатора


Код

def plot_boundary(clf, X, y, grid_step=.01, poly_featurizer=None): x_min, x_max = X[:, 0].min() - .1, X[:, 0].max() + .1 y_min, y_max = X[:, 1].min() - .1, X[:, 1].max() + .1 xx, yy = np.meshgrid(np.arange(x_min, x_max, grid_step), np.arange(y_min, y_max, grid_step)) # каждой точке в сетке x # ставим в соответствие свой цвет Z = clf.predict(poly_featurizer.transform(np.c_)) Z = Z.reshape(xx.shape) plt.contour(xx, yy, Z, cmap=plt.cm.Paired)


Полиномиальными признаками до степени для двух переменных и мы называем следующие:



Например, для это будут следующие признаки:



Нарисовав треугольник Пифагора, Вы сообразите, сколько таких признаков будет для и вообще для любого .
Попросту говоря, таких признаков экспоненциально много, и строить, скажем, для 100 признаков полиномиальные степени 10 может оказаться затратно (а более того, и не нужно).


Создадим объект sklearn, который добавит в матрицу полиномиальные признаки вплоть до степени 7 и обучим логистическую регрессию с параметром регуляризации . Изобразим разделяющую границу.
Также проверим долю правильных ответов классификатора на обучающей выборке. Видим, что регуляризация оказалась слишком сильной, и модель "недообучилась". Доля правильных ответов классификатора на обучающей выборке оказалась равной 0.627.


Код

poly = PolynomialFeatures(degree=7) X_poly = poly.fit_transform(X)
C = 1e-2 logit = LogisticRegression(C=C, n_jobs=-1, random_state=17) logit.fit(X_poly, y) plot_boundary(logit, X, y, grid_step=.01, poly_featurizer=poly) plt.scatter(X, X, c="green", label="Выпущен") plt.scatter(X, X, c="red", label="Бракован") plt.xlabel("Тест 1") plt.ylabel("Тест 2") plt.title("2 теста микрочипов. Логит с C=0.01") plt.legend(); print("Доля правильных ответов классификатора на обучающей выборке:", round(logit.score(X_poly, y), 3))


Увеличим до 1. Тем самым мы ослабляем регуляризацию, теперь в решении значения весов логистической регрессии могут оказаться больше (по модулю), чем в прошлом случае. Теперь доля правильных ответов классификатора на обучающей выборке – 0.831.


Код

C = 1 logit = LogisticRegression(C=C, n_jobs=-1, random_state=17) logit.fit(X_poly, y) plot_boundary(logit, X, y, grid_step=.005, poly_featurizer=poly) plt.scatter(X, X, c="green", label="Выпущен") plt.scatter(X, X, c="red", label="Бракован") plt.xlabel("Тест 1") plt.ylabel("Тест 2") plt.title("2 теста микрочипов. Логит с C=1") plt.legend(); print("Доля правильных ответов классификатора на обучающей выборке:", round(logit.score(X_poly, y), 3))


Еще увеличим – до 10 тысяч. Теперь регуляризации явно недостаточно, и мы наблюдаем переобучение. Можно заметить, что в прошлом случае (при =1 и "гладкой" границе) доля правильных ответов модели на обучающей выборке не намного ниже, чем в 3 случае, зато на новой выборке, можно себе представить, 2 модель сработает намного лучше.
Доля правильных ответов классификатора на обучающей выборке – 0.873.


Код

C = 1e4 logit = LogisticRegression(C=C, n_jobs=-1, random_state=17) logit.fit(X_poly, y) plot_boundary(logit, X, y, grid_step=.005, poly_featurizer=poly) plt.scatter(X, X, c="green", label="Выпущен") plt.scatter(X, X, c="red", label="Бракован") plt.xlabel("Тест 1") plt.ylabel("Тест 2") plt.title("2 теста микрочипов. Логит с C=10k") plt.legend(); print("Доля правильных ответов классификатора на обучающей выборке:", round(logit.score(X_poly, y), 3))


Чтоб обсудить результаты, перепишем формулу для функционала, который оптимизируется в логистической регрессии, в таком виде:





Промежуточные выводы :



Настройка параметра регуляризации


Теперь найдем оптимальное (в данном примере) значение параметра регуляризации . Сделать это можно с помощью LogisticRegressionCV – перебора параметров по сетке с последующей кросс-валидацией. Этот класс создан специально для логистической регрессии (для нее известны эффективные алгоритмы перебора параметров), для произвольной модели мы бы использовали GridSearchCV, RandomizedSearchCV или, например, специальные алгоритмы оптимизации гиперпараметров, реализованные в hyperopt.


Код

skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=17) c_values = np.logspace(-2, 3, 500) logit_searcher = LogisticRegressionCV(Cs=c_values, cv=skf, verbose=1, n_jobs=-1) logit_searcher.fit(X_poly, y)


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


Выделим участок с "лучшими" значениями C.


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

4. Где логистическая регрессия хороша и где не очень

Анализ отзывов IMDB к фильмам

Будем решать задачу бинарной классификации отзывов IMDB к фильмам. Имеется обучающая выборка с размеченными отзывами, по 12500 отзывов известно, что они хорошие, еще про 12500 – что они плохие. Здесь уже не так просто сразу приступить к машинному обучению, потому что готовой матрицы нет – ее надо приготовить. Будем использовать самый простой подход – мешок слов ("Bag of words"). При таком подходе признаками отзыва будут индикаторы наличия в нем каждого слова из всего корпуса, где корпус – это множество всех отзывов. Идея иллюстрируется картинкой


from __future__ import division, print_function # отключим всякие предупреждения Anaconda import warnings warnings.filterwarnings("ignore") import numpy as np import matplotlib.pyplot as plt %matplotlib inline import seaborn as sns import numpy as np from sklearn.datasets import load_files from sklearn.feature_extraction.text import CountVectorizer, TfidfTransformer, TfidfVectorizer from sklearn.linear_model import LogisticRegression from sklearn.svm import LinearSVC

Загрузим данные отсюда (краткое описание - ). В обучающей и тестовой выборках по 12500 тысяч хороших и плохих отзывов к фильмам.


reviews_train = load_files("YOUR PATH") text_train, y_train = reviews_train.data, reviews_train.target
print("Number of documents in training data: %d" % len(text_train)) print(np.bincount(y_train))
# поменяйте путь к файлу reviews_test = load_files("YOUR PATH") text_test, y_test = reviews_test.data, reviews_test.target print("Number of documents in test data: %d" % len(text_test)) print(np.bincount(y_test))

Пример плохого отзыва:


"Words can\"t describe how bad this movie is. I can\"t explain it by writing only. You have too see it for yourself to get at grip of how horrible a movie really can be. Not that I recommend you to do that. There are so many clich\xc3\xa9s, mistakes (and all other negative things you can imagine) here that will just make you cry. To start with the technical first, there are a LOT of mistakes regarding the airplane. I won\"t list them here, but just mention the coloring of the plane. They didn\"t even manage to show an airliner in the colors of a fictional airline, but instead used a 747 painted in the original Boeing livery. Very bad. The plot is stupid and has been done many times before, only much, much better. There are so many ridiculous moments here that i lost count of it really early. Also, I was on the bad guys\" side all the time in the movie, because the good guys were so stupid. "Executive Decision" should without a doubt be you\"re choice over this one, even the "Turbulence"-movies are better. In fact, every other movie in the world is better than this one."

Пример хорошего отзыва:


"Everyone plays their part pretty well in this "little nice movie". Belushi gets the chance to live part of his life differently, but ends up realizing that what he had was going to be just as good or maybe even better. The movie shows us that we ought to take advantage of the opportunities we have, not the ones we do not or cannot have. If U can get this movie on video for around $10, it\xc2\xb4d be an investment!"

Простой подсчет слов

Составим словарь всех слов с помощью CountVectorizer. Всего в выборке 74849 уникальных слов. Если посмотреть на примеры полученных "слов" (лучше их называть токенами), то можно увидеть, что многие важные этапы обработки текста мы тут пропустили (автоматическая обработка текстов – это могло бы быть темой отдельной серии статей).


Код

cv = CountVectorizer() cv.fit(text_train) print(len(cv.vocabulary_)) #74849
print(cv.get_feature_names()[:50]) print(cv.get_feature_names())

["00", "000", "0000000000001", "00001", "00015", "000s", "001", "003830", "006", "007", "0079", "0080", "0083", "0093638", "00am", "00pm", "00s", "01", "01pm", "02", "020410", "029", "03", "04", "041", "05", "050", "06", "06th", "07", "08", "087", "089", "08th", "09", "0f", "0ne", "0r", "0s", "10", "100", "1000", "1000000", "10000000000000", "1000lb", "1000s", "1001", "100b", "100k", "100m"]
["pincher", "pinchers", "pinches", "pinching", "pinchot", "pinciotti", "pine", "pineal", "pineapple", "pineapples", "pines", "pinet", "pinetrees", "pineyro", "pinfall", "pinfold", "ping", "pingo", "pinhead", "pinheads", "pinho", "pining", "pinjar", "pink", "pinkerton", "pinkett", "pinkie", "pinkins", "pinkish", "pinko", "pinks", "pinku", "pinkus", "pinky", "pinnacle", "pinnacles", "pinned", "pinning", "pinnings", "pinnochio", "pinnocioesque", "pino", "pinocchio", "pinochet", "pinochets", "pinoy", "pinpoint", "pinpoints", "pins", "pinsent"]


Закодируем предложения из текстов обучающей выборки индексами входящих слов. Используем разреженный формат. Преобразуем так же тестовую выборку.


X_train = cv.transform(text_train) X_test = cv.transform(text_test)

Обучим логистическую регрессию и посмотрим на доли правильных ответов на обучающей и тестовой выборках. Получается, на тестовой выборке мы правильно угадываем тональность примерно 86.7% отзывов.


Код

%%time logit = LogisticRegression(n_jobs=-1, random_state=7) logit.fit(X_train, y_train) print(round(logit.score(X_train, y_train), 3), round(logit.score(X_test, y_test), 3))


Коэффициенты модели можно красиво отобразить.


Код визуализации коэффициентов модели

def visualize_coefficients(classifier, feature_names, n_top_features=25): # get coefficients with large absolute values coef = classifier.coef_.ravel() positive_coefficients = np.argsort(coef)[-n_top_features:] negative_coefficients = np.argsort(coef)[:n_top_features] interesting_coefficients = np.hstack() # plot them plt.figure(figsize=(15, 5)) colors = ["red" if c < 0 else "blue" for c in coef] plt.bar(np.arange(2 * n_top_features), coef, color=colors) feature_names = np.array(feature_names) plt.xticks(np.arange(1, 1 + 2 * n_top_features), feature_names, rotation=60, ha="right");
def plot_grid_scores(grid, param_name): plt.plot(grid.param_grid, grid.cv_results_["mean_train_score"], color="green", label="train") plt.plot(grid.param_grid, grid.cv_results_["mean_test_score"], color="red", label="test") plt.legend();
visualize_coefficients(logit, cv.get_feature_names())




Подберем коэффициент регуляризации для логистической регрессии. Используем sklearn.pipeline, поскольку CountVectorizer правильно применять только на тех данных, на которых в текущий момент обучается модель (чтоб не "подсматривать" в тестовую выборку и не считать по ней частоты вхождения слов). В данном случае pipeline задает последовательность действий: применить CountVectorizer, затем обучить логистическую регрессию. Так мы поднимаем долю правильных ответов до 88.5% на кросс-валидации и 87.9% – на отложенной выборке.


Код

from sklearn.pipeline import make_pipeline text_pipe_logit = make_pipeline(CountVectorizer(), LogisticRegression(n_jobs=-1, random_state=7)) text_pipe_logit.fit(text_train, y_train) print(text_pipe_logit.score(text_test, y_test)) from sklearn.model_selection import GridSearchCV param_grid_logit = {"logisticregression__C": np.logspace(-5, 0, 6)} grid_logit = GridSearchCV(text_pipe_logit, param_grid_logit, cv=3, n_jobs=-1) grid_logit.fit(text_train, y_train) grid_logit.best_params_, grid_logit.best_score_ plot_grid_scores(grid_logit, "logisticregression__C") grid_logit.score(text_test, y_test)




Теперь то же самое, но со случайным лесом. Видим, что с логистической регрессией мы достигаем большей доли правильных ответов меньшими усилиями. Лес работает дольше, на отложенной выборке 85.5% правильных ответов.


Код для обучения случайного леса

from sklearn.ensemble import RandomForestClassifier forest = RandomForestClassifier(n_estimators=200, n_jobs=-1, random_state=17) forest.fit(X_train, y_train) print(round(forest.score(X_test, y_test), 3))

XOR-проблема

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


Линейные методы классификации строят все же очень простую разделяющую поверхность – гиперплоскость. Самый известный игрушечный пример, в котором классы нельзя без ошибок поделить гиперплоскостью (то есть прямой, если это 2D), получил имя "the XOR problem".


XOR – это "исключающее ИЛИ", булева функция со следующей таблицей истинности:



XOR дал имя простой задаче бинарной классификации, в которой классы представлены вытянутыми по диагоналям и пересекающимися облаками точек.


Код, рисующий следующие 3 картинки

# порождаем данные rng = np.random.RandomState(0) X = rng.randn(200, 2) y = np.logical_xor(X[:, 0] > 0, X[:, 1] > 0)
plt.scatter(X[:, 0], X[:, 1], s=30, c=y, cmap=plt.cm.Paired);
def plot_boundary(clf, X, y, plot_title): xx, yy = np.meshgrid(np.linspace(-3, 3, 50), np.linspace(-3, 3, 50)) clf.fit(X, y) # plot the decision function for each datapoint on the grid Z = clf.predict_proba(np.vstack((xx.ravel(), yy.ravel())).T)[:, 1] Z = Z.reshape(xx.shape) image = plt.imshow(Z, interpolation="nearest", extent=(xx.min(), xx.max(), yy.min(), yy.max()), aspect="auto", origin="lower", cmap=plt.cm.PuOr_r) contours = plt.contour(xx, yy, Z, levels=, linewidths=2, linetypes="--") plt.scatter(X[:, 0], X[:, 1], s=30, c=y, cmap=plt.cm.Paired) plt.xticks(()) plt.yticks(()) plt.xlabel(r"$$") plt.ylabel(r"$$") plt.axis([-3, 3, -3, 3]) plt.colorbar(image) plt.title(plot_title, fontsize=12);
plot_boundary(LogisticRegression(), X, y, "Logistic Regression, XOR problem")
from sklearn.preprocessing import PolynomialFeatures from sklearn.pipeline import Pipeline
logit_pipe = Pipeline([("poly", PolynomialFeatures(degree=2)), ("logit", LogisticRegression())])
plot_boundary(logit_pipe, X, y, "Logistic Regression + quadratic features. XOR problem")


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


А вот если на вход подать полиномиальные признаки, в данном случае до 2 степени, то проблема решается.


Здесь логистическая регрессия все равно строила гиперплоскость, но в 6-мерном пространстве признаков и . В проекции на исходное пространство признаков граница получилась нелинейной.


На практике полиномиальные признаки действительно помогают, но строить их явно – вычислительно неэффективно. Гораздо быстрее работает SVM с ядровым трюком. При таком подходе в пространстве высокой размерности считается только расстояние между объектами (задаваемое функцией-ядром), а явно плодить комбинаторно большое число признаков не приходится. Про это подробно можно почитать в курсе Евгения Соколова (математика уже серьезная).

5. Кривые валидации и обучения

Мы уже получили представление о проверке модели, кросс-валидации и регуляризации.
Теперь рассмотрим главный вопрос:


Если качество модели нас не устраивает, что делать?

  • Сделать модель сложнее или упростить?
  • Добавить больше признаков?
  • Или нам просто нужно больше данных для обучения?

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


Будем работать со знакомыми данными по оттоку клиентов телеком-оператора.


Импорт библиотек и чтение данных

from __future__ import division, print_function # отключим всякие предупреждения Anaconda import warnings warnings.filterwarnings("ignore") %matplotlib inline from matplotlib import pyplot as plt import seaborn as sns import numpy as np import pandas as pd from sklearn.preprocessing import PolynomialFeatures from sklearn.pipeline import Pipeline from sklearn.preprocessing import StandardScaler from sklearn.linear_model import LogisticRegression, LogisticRegressionCV, SGDClassifier from sklearn.model_selection import validation_curve data = pd.read_csv("../../data/telecom_churn.csv").drop("State", axis=1) data["International plan"] = data["International plan"].map({"Yes": 1, "No": 0}) data["Voice mail plan"] = data["Voice mail plan"].map({"Yes": 1, "No": 0}) y = data["Churn"].astype("int").values X = data.drop("Churn", axis=1).values


Логистическую регрессию будем обучать стохастическим градиентным спуском. Пока объясним это тем, что так быстрее, но далее в программе у нас отдельная статья про это дело. Построим валидационные кривые, показывающие, как качество (ROC AUC) на обучающей и проверочной выборке меняется с изменением параметра регуляризации.


Код

alphas = np.logspace(-2, 0, 20) sgd_logit = SGDClassifier(loss="log", n_jobs=-1, random_state=17) logit_pipe = Pipeline([("scaler", StandardScaler()), ("poly", PolynomialFeatures(degree=2)), ("sgd_logit", sgd_logit)]) val_train, val_test = validation_curve(logit_pipe, X, y, "sgd_logit__alpha", alphas, cv=5, scoring="roc_auc") def plot_with_err(x, data, **kwargs): mu, std = data.mean(1), data.std(1) lines = plt.plot(x, mu, "-", **kwargs) plt.fill_between(x, mu - std, mu + std, edgecolor="none", facecolor=lines.get_color(), alpha=0.2) plot_with_err(alphas, val_train, label="training scores") plot_with_err(alphas, val_test, label="validation scores") plt.xlabel(r"$\alpha$"); plt.ylabel("ROC AUC") plt.legend();


Тенденция видна сразу, и она очень часто встречается.

    Для простых моделей тренировочная и валидационная ошибка находятся где-то рядом, и они велики. Это говорит о том, что модель недообучилась : то есть она не имеет достаточное кол-во параметров.

    Для сильно усложненных моделей тренировочная и валидационная ошибки значительно отличаются. Это можно объяснить переобучением : когда параметров слишком много либо не хватает регуляризации, алгоритм может "отвлекаться" на шум в данных и упускать основной тренд.

Сколько нужно данных?


Известно, что чем больше данных использует модель, тем лучше. Но как нам понять в конкретной ситуации, помогут ли новые данные? Скажем, целесообразно ли нам потратить N\$ на труд асессоров, чтобы увеличить выборку вдвое?


Поскольку новых данных пока может и не быть, разумно поварьировать размер имеющейся обучающей выборки и посмотреть, как качество решения задачи зависит от объема данных, на котором мы обучали модель. Так получаются кривые обучения (learning curves ).


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


Давайте посмотрим, что мы получим для линейной модели. Коэффициент регуляризации выставим большим.


Код

from sklearn.model_selection import learning_curve def plot_learning_curve(degree=2, alpha=0.01): train_sizes = np.linspace(0.05, 1, 20) logit_pipe = Pipeline([("scaler", StandardScaler()), ("poly", PolynomialFeatures(degree=degree)), ("sgd_logit", SGDClassifier(n_jobs=-1, random_state=17, alpha=alpha))]) N_train, val_train, val_test = learning_curve(logit_pipe, X, y, train_sizes=train_sizes, cv=5, scoring="roc_auc") plot_with_err(N_train, val_train, label="training scores") plot_with_err(N_train, val_test, label="validation scores") plt.xlabel("Training Set Size"); plt.ylabel("AUC") plt.legend() plot_learning_curve(degree=2, alpha=10)


Типичная ситуация: для небольшого объема данных ошибки на обучающей выборке и в процессе кросс-валидации довольно сильно отличаются, что указывает на переобучение. Для той же модели, но с большим объемом данных ошибки "сходятся", что указывается на недообучение.


Если добавить еще данные, ошибка на обучающей выборке не будет расти, но с другой стороны, ошибка на тестовых данных не будет уменьшаться.


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


Что будет, если изменить коэффициент регуляризации (уменьшить до 0.05)?


Видим хорошую тенденцию – кривые постепенно сходятся, и если дальше двигаться направо (добавлять в модель данные), можно еще повысить качество на валидации.


А если усложнить модель ещё больше ()?


Проявляется переобучение – AUC падает как на обучении, так и на валидации.


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


Выводы по кривым валидации и обучения

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

6. Плюсы и минусы линейных моделей в задачах машинного обучения



  • Плохо работают в задачах, в которых зависимость ответов от признаков сложная, нелинейная
  • На практике предположения теоремы Маркова-Гаусса почти никогда не выполняются, поэтому чаще линейные методы работают хуже, чем, например, SVM и ансамбли (по качеству решения задачи классификации/регрессии)

7. Домашнее задание № 4

Актуальные домашние задания объявляются во время очередной сессии курса, следить можно в