Consideration of estimation of logistic regression as an optimization problem.pdf Логистическая регрессия является одним из распространенных методов классификации данных в разных областях [1. С. 181-242; 2. С. 283-369; 3. С. 152-214; 4. С. 321-329]. Основная цель логистической регрессии состоит в разделении множества исходных значений линейной границей (разделяющей прямой, плоскостью или гиперплоскостью) на две соответствующих двум заданным классам области. Логистическая регрессия прогнозирует вероятность некоторого события, находящуюся в пределах от 0 до 1 [4. С. 322-324]. Следует указать, что число классов L может быть больше двух, в этом случае они имеют мультиномиальную регрессию, которую можно построить с помощью L - 1 независимых логистических регрессий. Поэтому будем рассматривать задачу классификации с двумя классами. В настоящее время построение логистической регрессии обычно осуществляют в виде оптимизационной задачи, реализующей критерий максимального правдоподобия [5. С. 124]. Ряд авторов [6. С. 31-32; 7] предлагают оценивать коэффициенты разделяющего уравнения методом Ньютона-Рафсона. Однако исследования на модельных данных этого и других алгоритмов оценивания логистической регрессии, реализованных в статистических пакетах STATISTICA и SPSS, показали неустойчивость оценок. Этот недостаток отмечается также в [7]. Цель статьи - исследование особенностей процедуры оценивания логистической регрессии как экстремальной задачи, а также разработка устойчивого алгоритма вычисления коэффициентов разделяющего уравнения. ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА 2017 № 40 Управление, вычислительная техника и информатика Пусть имется выборка прецедентов (обучающая выборка): (xi,y), i = 1,2, ■..,n , (1) Л ( *Л ^ i1 x1 T Х 2 = T , X „ 1 x2 V 1 Xn2 М2 1m x X i 2 где xi = V Xim у X nm x, f 1 Л X i 2 - вектор значений i-го объекта, Vi xi1 = 1; X = V Xim у yi e {- 1; 1} - бинарная переменная, указывающая на принадлежность i-го объекта соответствующему классу, например, для определенности первому классу при yi = - 1 и второму - при yi = 1; M = m - 1 -число признаков у каждого объекта; n - количество наблюдений. Классификацию осуществляют с помощью функции h(x) =-1---, принимающей значе- 1 + exp{-b1 x} ния в интервале (0; 1). Пороговым значением является h(x) = 0,5 . В [6. С. 29-32] показано, что максимизация логарифма правдоподобия эквивалентна минимизации функционала Q(b) = > Ц1 + e-yibTx ) ^ min , (2) £i beRm где b T = (b1 b2 ... bm) - искомый вектор коэффициентов разделяющей линейной границы, описываемой в общем случае уравнением гиперплоскости m b1 + 1 bjXj = 0. j=2 Аналитически задача минимизации (2) не решается. Поэтому для оценивания вектора коэффициентов b используют итерационные алгоритмы спуска для решения экстремальных задач. В [6. С. 32-34] описан итерационный алгоритм Ньютона-Рафсона решения задачи (2). Он состоит в следующем. В качестве нулевого приближения можно взять решение задачи классификации методом многомерной линейной регрессии b(0) = (XT X)-1(XT Y). (3) Затем начинается итерационный процесс, на k-м шаге которого уточняется вектор коэффициентов b(k): b(k) = b(k-1) -hk(Q"(b(k-1)))-1 Q'(b(k-1)), где Q'(b(k)) - вектор первых производных (градиент) функционала Q(b) в точке b(k); Q''(b(k)) - матрица вторых производных (гессиан) функционала Q(b) в точке b(k); hk - величина шага, который можно положить равным 1, но его подбор на каждом шаге способен увеличить скорость сходимости. Метод Ньютона-Рафсона описан в ряде учебников по методам оптимизации, например в [8. С. 223-226]. Статистическое моделирование методом Монте-Карло [9] показало неустойчивость в некоторых случаях работы этого алгоритма, а также других алгоритмов, реализованных в статистических пакетах. Отметим основные причины этого. В итерационных алгоритмах поиска минимума требуется задать начальное приближение вектора b(0). Плохое начальное приближение, значительно отличающееся от точки локального минимума, может значительно ухудшить условия сходимости итерационного алгоритма. В результате статистического моделирования выяснилось, что использование в качестве начального приближения b(0) решения задачи классификации в виде линейной регрессии часто дает неудовлетворительный результат. Покажем это на примере. f 1 ^ fx ^ 11 Пример 1. Сгенерируем выборку прецедентов (x i, yi), x i = i = 1, 2,..., n, следую xi2 xi2 V xi3 J щим образом. Объем выборки n = 100, число признаков M = 2, первая половина выборки (i = 1,2,...,50) относится к первому классу (yi =-1), а вторая - ко второму классу (yi = 1). Пусть в каждом классе признаки являются взаимно независимыми гауссовскими случайными величинами. Это означает, что первый и второй классы представляют собой случайные выборки из двумерных гауссовских случайных векторов X2 = (XX3^) и X(2) = (X22), X32)) соответственно. Зададим у всех компонент одно и то же среднее квадратическое отклонение сX = 1,5 и математические ожидания: M [X21)] = 10, M[X3(1)] = 16 , M[X 22)] = 4 , M[X3(2)] = 4 . На рисунке 1 приведен результат классификации с помощью метода наименьших квадратов на основе (3). Разделяющая прямая не только не классифицирует выборку прецедентов на два класса, но и значительно удалена от всех объектов. Хотя классы в этом примере легкоразделимы. Такое неудачное начальное приближение значительно затрудняет работу итерационных алгоритмов решения задачи (2). V xi3 J Вторая причина неустойчивости решений задачи (2) состоит в том, что итерационные алгоритмы спуска рассчитаны на поиск локального минимума. А задача (2) при корректной классификации всех прецедентов имеет нижнюю грань на бесконечности, равную нулю. Причем качество построения разделяющей границы при формальном спуске уже не будет связано со значением целевой функции Q(b). Поэтому требуется корректировка стратегии спуска. Третья причина состоит в том, что при достаточно малых значениях целевой функции Q(b) некоторые величины - y i bT x i могут принимать очень большие значения, что приводит к неустранимым вычислительным погрешностям и даже к останову работы алгоритма из-за переполнения памяти. -80 х3 20 10 10 -20 -30 -40 50 -60 -70 •• * • * ... X* * V* * * 10 12 14 х2 Рис. 1. Пример построения разделяющей границы с помощью линейной регрессии Таким образом, приведенные недостатки непосредственного использования итерационных алгоритмов спуска требуют исследования особенностей задачи (2). 2. Исследование целевой функции задачи оценивания коэффициентов разделяющей границы по критерию максимального правдоподобия 1 + expj- y, £ bk Исследуем свойства целевой функции Q(b). Рассмотрим вначале произвольное i-е слагаемое целевой функции Q(b), т.е. функцию Z (b) = ln(1 + e-yib Xi) = ln к=1 Найдем ее частные производные: ybT x dz, y,x,ke к = 1, 2, ..., m . 1 + edbk dz. Так как все частные производные непрерывно дифференцируемы и -- Ф 0 , то функция (b) не имеет db j стационарных точек, а значит, у нее нет экстремумов. d2 z. . Они равны Далее найдем вторые производные dbkdbj ■ = yf xikxi,Ui(b), k, i =1,2, - •, m, dbkdbj y.b xi e где Ux (b) = - (1 + e-yb xi )2 Поскольку все вторые производные всюду непрерывны, то ■> 0. d2 Zi d2 Zi . Исследуем функцию dbkdbj dbjdbk zt (b) на выпуклость. Ее гессиан равен H (Zi) = (qj (b)), qj (b) = d2 Zi dbk dbj = yf xikxjui(bK k, j =1,2, •••, m. Воспользуемся критерием Сильвес тра и определим главные миноры Ax, A2, •.., Am гессиаn n на H(Q). Поскольку xn _ то A\ _ f xf\Ui (b) _ f U,. (b) > 0 как сумма строго положительных чисел. i_\ _ Остальные главные миноры с учетом выпуклости Q(b) равны fU, (b) fx^U, (b) • f^x^U, (b) i _\ i_\ i_\ f x,2U, (b) f x]2Ul (b) f xi2x,kU, (b) i_\ f xiU (b) Ak _ > 0, k _ 2, 3, m . i _\ i _\ f x,kU, (b) f x,kx,2U, (b) i_\ Рассмотрим частные случаи. Пусть m _ 2, n _ 2 . После преобразований получим A2 _ (x\2 - x22)2 U\ (b)U2 (b) > 0 . n-\ n Пусть m _ 2, n _ 3 . После преобразований получим A2 _ f f (xi 2 - x, 2 )2 U, (b)U, (b) > 0. i_\ j_i+\ Таким образом, квадратичная форма в случае нескольких наблюдений может быть положительно определенной или неотрицательно определенной. Это означает, что функция Q(b): " т» m - является выпуклой на К и не имеет максимумов; - может иметь стационарные точки, в которых могут быть локальные минимумы. Поскольку равенство нулю главных миноров будет только в некоторых точках, например в вырожденных случаях равенства всех измерений, то функция Q(b) может оказаться в окрестности стационарной точки строго выпуклой. Это означает, что в такой стационарной точке у функции Q(b) будет локальный минимум. Из этого можно сделать выводы. Использование численных методов спуска первого или второго порядков при наличии у целевой функции стационарных точек может привести к неправильному построению логистической регрессии по двум причинам. Во-первых, можно в качестве решения получить локальный минимум, который будет грубой оценкой коэффициентов уравнения логистической регрессии. Во-вторых, спуск может стремиться к нижней грани целевой функции, а вектор коэффициентов логистической регрессии будет оставаться грубым решением. 3. Устойчивый алгоритм вычисления коэффициентов разделяющего уравнения логистической регрессии i _\ i _\ f x,kU, (b) f x,kx,2U, (b) i_\ Рассмотрим частные случаи. Пусть m _ 2, n _ 2 . После преобразований получим A2 _ (x\2 - x22)2 U\ (b)U2 (b) > 0 . n-\ n Пусть m _ 2, n _ 3 . После преобразований получим A2 _ f f (xi 2 - x, 2 )2 U, (b)U, (b) > 0. i_\ j_i+\ Таким образом, квадратичная форма в случае нескольких наблюдений может быть положительно определенной или неотрицательно определенной. Это означает, что функция Q(b): " т» m - является выпуклой на К и не имеет максимумов; - может иметь стационарные точки, в которых могут быть локальные минимумы. Поскольку равенство нулю главных миноров будет только в некоторых точках, например в вырожденных случаях равенства всех измерений, то функция Q(b) может оказаться в окрестности стационарной точки строго выпуклой. Это означает, что в такой стационарной точке у функции Q(b) будет локальный минимум. Из этого можно сделать выводы. Использование численных методов спуска первого или второго порядков при наличии у целевой функции стационарных точек может привести к неправильному построению логистической регрессии по двум причинам. Во-первых, можно в качестве решения получить локальный минимум, который будет грубой оценкой коэффициентов уравнения логистической регрессии. Во-вторых, спуск может стремиться к нижней грани целевой функции, а вектор коэффициентов логистической регрессии будет оставаться грубым решением. 3. Устойчивый алгоритм вычисления коэффициентов разделяющего уравнения логистической регрессии Имеем выборку прецедентов (\), в которой первые N\ объектов относятся к первому классу (у, _ - \), а остальные - ко второму классу (y, _ \). С учетом вышеизложенного приведем описание алгоритма. Он состоит в следующем. L Найдем середину x(0) отрезка, соединяющего центры тяжести первого и второго классов i _\ ! xi , x 3 _ _-f x, , N\ t! f x n - N\ ,_N\ +\ 3. Через точку x(0) строим гиперплоскость S: s1 + s2x2 +... + smxm = 0, ортогональную вектору d = xl/J _ Xw. 4. Получаем начальную оценку b(0) вектора разделяющей гиперплоскости как b(0) = s /| s|, где s -вектор коэффициентов гиперплоскости S. При этом длина b(0) = 1. Зададим параметр р 0 = 1. 5. С помощью итерационного алгоритма ищем с заданной точностью в решение задачи (2). На любом k-м шаге алгоритма, многократно генерируя случайный вектор p(l), решаем задачу '22 = X (2) _ X (j) Q(b(k,t)) ^ mmm, b(k >CR m b(kJ) = a(b(kl _j) + hk p(l)) k (k,1) = b(k_1). (k ,i) = P k _1, b(k ,j) = b Если требуемая точность решения будет достигнута (Q(b(k,l)) 1 - коэффициент растяжения и возвращаемся к п. 5 (k := k +1). Следует отметить, что слишком малое значение параметра точности в может привести к большим длинам оценок векторов b(k,l) и, как указано выше, к потере точности алгоритма. 4. Вычислительные эксперименты Проведем с помощью метода статистических испытаний Монте-Карло несколько экспериментов для исследования эффективности предложенного алгоритма по сравнению с алгоритмом «Quasi-Newton», реализованным в статистическом пакете STATISTICA. Он показал наилучшие результаты оценивания среди алгоритмов логистической регрессии в STATISTICA. Число опытов L = 200 . Выборка прецедентов имеет тот же вид, что и в примере 1, различия заключаются в объемах выборок первого N1 и второго N2 классов (N1 + N2 = n) и математических ожиданиях M[X21) ], M[X3(1) ], M[X22) ], M[X3(2) ]. Поскольку с точностью до постоянного множителя векторов коэффициентов b может быть бесконечное множество, то для удобства сравнения результаты оценивания будем задавать в нормированном виде с единичной длиной. Пример 2. Рассмотрим случай одинаковых объемов выборок в классах. Зададим N1 = N2 = 20. Результаты оценивания коэффициентов логистической регрессии по предложенному алгоритму и алгоритму «Quasi-Newton» в пакете STATISTICA приведены в табл. 1. Т а б л и ц а 1 Результаты оценивания логистической регрессии по предложенному алгоритму и алгоритму «Quasi-Newton» в пакете STATISTICA для случая N = N2 = 20 Матожидания признаков Параметр Предложенн^1й алгоритм Алгоритм «Quasi-Newton» b1 b2 b3 b1 b2 b3 M[ X f] = 6, M[ X 3(1)] = 3, M[ x 22)] = _6, M[ X 3(2)] = _3 Pj 0,00000 -0,89443 -0,44721 0,00000 -0,89443 -0,44721 bj -0,00925 -0,88467 -0,46612 -0,13761 -0,78647 -0,28092 s(bj) 0,26070 0,10410 0,13367 0,47889 0,21434 0,29080 Abj -0,00925 0,00975 -0,01891 -0,13761 0,10795 0,16629 M[ X f] = 16, M[ X 3(1)] = 13, M[ X 22)] = 4, M[ X 3(2)] = 7 Pj 0,99723 -0,06648 -0,03324 0,99723 -0,06648 -0,03324 bj 0,99721 -0,06691 -0,03297 0,99701 -0,04299 -0,05718 s(bj) 0,00047 0,00958 0,00923 0,00063 0,03115 0,03347 Abj -0,00002 -0,00042 0,00027 -0,00023 0,02350 -0,02394 Матожидания Параметр Предложенный алгоритм Алгоритм «Quasi-Newton» признаков b1 b2 b3 b1 b2 b3 M[ X «] = 26, Pj 0,99931 -0,03331 -0,01666 0,99931 -0,03331 -0,01666 M[X3a)] = 23 , bj 0,99930 -0,03345 -0,01648 0,99917 -0,03187 -0,01801 M[ X 22)] = 14, s(bj) 0,00009 0,00428 0,00429 0,00060 0,01393 0,01250 M[ X 3(2)] = 17 Abj 0,00000 -0,00014 0,00018 -0,00014 0,00144 -0,00135 Пример 3. Рассмотрим случай разных объемов выборок в классах. Зададим N1 = 20, N2 = 40 . Результаты оценивания коэффициентов логистической регрессии по предложенному алгоритму и алгоритму «Quasi-Newton» в пакете STATISTICA приведены в табл. 2. Т а б л и ц а 2 Результаты оценивания логистической регрессии по предложенному алгоритму и алгоритму «Quasi-Newton» в пакете STATISTICA для случая N1 = 20, N2 = 40 Матожидания признаков Параметр Предложенный алгоритм Алгоритм «Quasi-Newton» b1 b2 b3 b1 b2 b3 M[ X f] = 6, M[ X f] = 3, M[ x 22)] = -6, M[ X 3(2)] = -3 Pj 0,00000 -0,89443 -0,44721 0,00000 -0,89443 -0,44721 bj 0,09740 -0,89644 -0,43234 -0,10119 -0,73443 -0,18157 s(bj) 0,25587 0,08658 0,14076 0,55009 0,27014 0,41928 Abj 0,09740 -0,00201 0,01488 -0,10119 0,15999 0,26565 M[ X f] = 16, M[ X 3(1)] = 13, M[ X 22)] = 4, M[ X 3(2)] = 7 Pj 0,99723 -0,06648 -0,03324 0,99723 -0,06648 -0,03324 bj 0,99724 -0,06699 -0,03194 0,99702 -0,03436 -0,06328 s(bj) 0,00035 0,00852 0,00881 0,00068 0,03844 0,03698 Abj 0,00001 -0,00051 0,00130 -0,00022 0,03213 -0,03004 M[ X f] = 26, M[X3a)] = 23 , M[ X 22)] = 14, M[ X 3(2)] = 17 Pj 0,99931 -0,03331 -0,01666 0,99931 -0,03331 -0,01666 bj 0,99931 -0,03303 -0,01676 0,99894 -0,02170 -0,02813 s(bj) 0,00010 0,00499 0,00493 0,00090 0,02441 0,02436 Abj 0,00001 0,00028 -0,00011 -0,00037 0,01161 -0,01148 В таблицах 1 и 2 приняты следующие обозначения: - теоретические значения коэффициентов логистической регрессии; bj - средние значения оценок коэффициентов логистической регрессии; s(bj) - средние квадратические значения ошибок оценок коэффициентов логистической регрессии; Abj -средние значения ошибок оценок коэффициентов логистической регрессии. Анализ таблиц показывает, что предложенный алгоритм значительно выигрывает по сравнению с известным. Заключение Исследования показали вычислительную неустойчивость алгоритмов оценивания логистической регрессии, реализующих критерий максимального правдоподобия. С целью повышения точности оценок предложен новый итерационный алгоритм. Он включает последовательность задач минимизации нулевого порядка на основе случайного поиска. На каждом k-м шаге длина вектора коэффициентов b(k) является фиксированной. Ее постепенно увеличивают до достижения целевой функцией требуемого значения. Результаты статистического моделирования показали более высокую точность оценок логистической регрессии предложенного алгоритма по сравнению с известными.

Скачать электронную версию публикации