Поиск параметрической регрессионной модели в индуктивно заданном множестве Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 12, № 1, 2007

ПОИСК ПАРАМЕТРИЧЕСКОЙ РЕГРЕССИОННОЙ МОДЕЛИ В ИНДУКТИВНО ЗАДАННОМ

МНОЖЕСТВЕ*

В. В. Стрижов

Вычислительный центр имени А. А. Дородницына РАН, Москва, Россия

e-mail: strijov@ccas.ru

A procedure of the search for a parametric regression model in a model set is described. The model set is a set of superpositions of the given smooth functions. The models' parameters density estimates are used for the search. To illustrate applicability of the approach the problem of the pressure variation in a spray chamber of the combustion engine is examined.

Введение

Проблема отыскания оптимальной параметрической регрессионной модели имеет большую историю, однако продолжает оставаться одной их самых актуальных в области распознавания образов. А.Г. Ивахненко, еще в 1968 году, предложил метод группового учета аргументов [1]. Согласно этому методу модель, доставляющая наилучшее приближение, отыскивается во множестве последовательно порождаемых моделей. В частности, для построения моделей как суперпозиций функций использовались полиномиальные функции, ряды Фурье и некоторые другие функции. А.Г. Ивахненко и его ученики создали ряд алгоритмов синтеза моделей и предложили методы оценки качества моделей.

При порождении конкурирующих моделей появляется задача определения значимости элементов модели. В работе К. Бишопа [2] предложен метод анализа распределения параметров однослойных нейронных сетей посредством гиперпараметров, т. е. параметров аппроксимирующих функций. Для каждого элемента сети оценивается плотность гауссов-ского распределения его параметров и делается вывод о том, насколько информативен данный элемент исследуемой регрессионной модели.

Ле Кюн предложил метод для модификации моделей, называемый методом оптимального отсечения (optimal brain damage) [3]. Этот метод состоит в исключении некоторых, наименее информативных, элементов регрессионной модели с тем условием, что при этом качество аппроксимации уменьшается незначительно. При исключении отдельных элементов модели становится возможным оценить их вклад по значениям заданной функции качества аппроксимации.

Проблема сравнения и выбора регрессионных моделей получила новое развитие после ряда публикаций Д. МакКая [4-6], предложившего при выборе модели из заданного мно-

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 07-07-00181).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

жества использовать не информационные критерии, например AIC — Akakie Information Criterion, а двухуровневый байесовский вывод и правило Оккама. На первом уровне вывода вычисляются плотности вероятностей распределения параметров каждой модели из заданного множества, на втором — правдоподобие моделей. Правило Оккама состоит в том, что вероятность выбора более сложной модели меньше, чем более простой, при сравнимом значении функции качества аппроксимации.

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

Поставим задачу нахождения регрессионной модели нескольких свободных переменных следующим образом. Задана выборка — множество {xi,...,Xn|x £ RM} значений свободных переменных и множество {у1,...,у^|у £ R} соответствующих им значений зависимой переменной. Обозначим оба эти множества как множество исходных данных D. Также задано множество G = {д|д : R х ... х R —> R} гладких параметрических функций д = g(b, •, •,..., •). Первый аргумент функции д — вектор-строка параметров b, последующие — переменные из множества действительных чисел, рассматриваемые как элементы вектора свободных переменных. Рассмотрим произвольную суперпозицию, состоящую из не более чем r функций д. Эта суперпозиция задает параметрическую регрессионную модель f = f (w, x). Модель f зависит от вектора свободных переменных x и вектора параметров w. Вектор w £ состоит из присоединенных вектор-параметров функций д1,...,дг, т.е. w = b1.b2.....br, где . — знак присоединения векторов. Обозначим через Ф = {fi} множество всех суперпозиций, индуктивно порожденное элементами множества G.

Требуется найти модель fi, которая доставляет максимум функционалаp(w|D, а, в, fi). Этот функционал, определяемый далее, включает искомую модель fi(w, x) и ее дополнительные параметры а и в.

1. Выбор регрессионных моделей и гипотеза порождения данных

Общий подход к сравнению нелинейных моделей заключается в следующем. Рассмотрим набор конкурирующих моделей /1,...,/м• Априорная вероятность модели /г определена как Р(/г). При появлении данных О апостериорная вероятность модели Р(/г|О) может быть найдена по теореме Байеса

P (fi|D)

P (fi )p(D|fi)

м :

Ep(D|fj )P (fj) j=i

где p(D| fi) — функция соответствия модели данным. Знаменатель дроби обеспечивает

M

выполнение условия У] P(fi|D) = 1.

i= 1

Вероятности моделей f1 и f2, параметры которых идентифицированы по данным D, сравнимы как

P(fi|D) P(fi)p(D|fi)

P (f2|D) P (f2)p(D|f2)'

(1)

Отношение р(О|/1)/р(О|/) есть отношение правдоподобия моделей. Отношение Р(/1 )/Р(/2) является априорной оценкой предпочтения одной модели другой. При моделировании отдается предпочтение наиболее простым и устойчивым моделям. Если априорные оценки Р(/¿) моделей одинаковы, т. е. нет причины предпочитать одну модель другой, то их необходимо сравнивать по значениям р(О|/):

Апостериорная плотность распределения параметров w функции / при заданной выборке О равна

рМД/,) = , (2)

где р^|/) — априорно заданная плотность вероятности параметров начального приближения; р(О^,/г) — функция правдоподобия параметров модели, а знаменатель р(О|/) обеспечивает выполнение условия /p(w|D,/¿)dw = 1. Он задан интегралом в пространстве параметров /р^'|О, Формулы (2) и (1) называются формулами байесовского вывода первого и второго уровня.

Рассмотрим регрессию у = /¿(Ь, х) + V с аддитивным гауссовским шумом с дисперсией аи и нулевым матожиданием. Тогда плотность вероятности появления данных опреде-

' V

ляется как

/I , лд /П| n n exp(-eED(D|w,fi))

Zd (в)

где в = l/^. Нормирующий множитель Zd (в) задан выражением

■ \ N

Zd (в) = ( j) 2 , (3)

а взвешенный функционал ошибки в пространстве данных

в^ = 2 £(/.(х») - У™)2- (4)

П=1

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

ехр(-аЕ^ М/Ш РМа^ = -

Zw (а)

где а — обратная дисперсия распределения параметров, а = а нормирующая константа Zw зависит от дисперсии распределения параметров как

ж

Zw(а) = 0т) 2 • (5)

Требование к малым значениям параметров [7] предполагает гауссовское априорное распределение с нулевым средним:

р(') = zWехР (-1!М|2)•

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

а

аEw = - |М|2. (6)

Эта ошибка регуляризирует параметры, начисляя штраф за их чрезмерно большие значения.

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

рНД^) =-—---.

р(Ща,в)

Записывая функцию ошибки в виде Б(') = аEw + вЕд, получаем

( ш а г\ ехР(-Б ('|/г)) РЫО,а,в,/г) = -—(—^-, (7)

Zs(а, в)

где Zs — нормирующий множитель.

2. Нахождение параметров модели

Рассмотрим итеративный алгоритм для определения оптимальных параметров ' и гиперпараметров а, в при заданной модели /г. Корректный подход заключается в интегрировании всех неизвестных параметров и гиперпараметров. Апостериорное распределение параметров определяется как

р('Щ = Цр(',а,в|D)dаdв = Цр(^а, в^)р(а, в|D)dаdв, (8)

что требует интегрирования апостериорного распределения параметров p(w|а, в, D) по пространству, размерность которого равна количеству параметров. Вычислительная сложность этого интегрирования весьма велика. Интеграл может быть упрощен при подходящем выборе начальных значений гиперпараметров.

Приближение интеграла заключается в том, что апостериорная плотность распределения гиперпараметров р(а, вимеет выраженный пик в окрестности наиболее правдоподобных значений гиперпараметров амр, вМР • Это приближение известно как аппроксимация Лапласа [8]. При таком допущении интеграл (8) упрощается до

рМЯ) « рНаМР,вМР,£)^Iр(а, в|D)dadв ~ p(w|амр, вмр, D).

Необходимо найти значения гиперпараметров, которые оптимизируют апостериорную плотность вероятности параметров, а затем выполнить все остальные расчеты, включающие р^|Д) при фиксированных значениях гиперпараметров.

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

БЫ « Б(wмр) + - wMP)TА^ - wмр). (9)

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

дБ^мр) _ 0 ди^

для всех значений Матрица А — это матрица Гессе функции ошибок:

А _ У2Б^мр) _ вV2ED(wмр) + а!.

Обозначим первое слагаемое правой части через Н, тогда А _ Н + а!. Подставив полученное приближенное значение Бв (7) и обозначив Дw _ w — wмр, получим

р^а,в^) _ -^ехр ( —Б^мр) — 2ДwTАД^ .

V 2 /

Оценим нормирующую константу "б, необходимую для аппроксимации кривой Гаусса, как

¿б _ ехр(—Б^мр))(2п)^(ае1 А)-2. (10)

Максимизируем функцию р(Д|а,в), изменяя значения гиперпараметров а и в• Это можно выполнить, интегрируя функцию плотности вероятности данных по пространству параметров w:

р(П^,в а, в )рМа,в )dw р(Б^,а,в )p(w|а)dw, (11)

где второй интеграл справедлив по причине того, что распределение параметров не зависит от дисперсии шума в силу гипотезы о гауссовском распределении шума. Для упрощения вычислений мы допускаем, что распределение р(а,в) является равномерным. Используя (4), (6), запишем (11) в виде

1 1 /* р(Д|а,в) _ ехр(—Б^^

Я в (в) Я в (а) }

Из (3), (5), (10) и предыдущего выражения получим

W

N

N.

1пр(£|а, в) = -а£$Р - в^Г - " 1п IА| + — 1п а + — 1пв--1п (2п).

2

2

2

2

(12)

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

1на| = Ла МП л +

^ 1 а/ ^ 4^ = 1 / ,=1 ,= 1

В этом выражении Л1;..., Лш — собственные значения матрицы Н. Так как функция ошибки на данных не является квадратичной функцией параметров, как при линейной или ИВЕ регрессии, непосредственно оптимизировать величину а невозможно, гессиан Н не является константой, а зависит от параметров Поскольку мы принимаем А = Н + а/ для вектора 'мр, который зависит от выбора а, собственные значения Н косвенным образом зависят от а. Таким образом, формула (13) игнорирует параметры модели. С использованием этого приближения, производная (12) с учетом а равна

Л ш Ш 1 * ^ 1п(Л, + а) = £ = 1г(А-1). (13)

Л о + а

1пР(^|а,в) = -ЕШР - 2Ё 1

W

2 ^ Л, + а 2а'

3 = 1

Приравнивая последнее выражение к нулю и преобразовывая его, получаем выражение

для а:

2аЕМР

W-

а

3=1 Лз + а'

(14)

3=1

Обозначим вычитаемое правой части через 7:

а

3 = 1

=1 Лз + а

Те компоненты суммы, в которых Л, ^ а, привносят вклад, близкий к единице, а те компоненты суммы, в которых 0 < Л, ^ а, привносят вклад, близкий к нулю. Таким образом 7 может быть интерпретирована как мера числа хорошо обусловленных параметров модели.

Для нахождения гиперпараметра в рассмотрим задачу оптимизации (12). Обозначим через собственное значение матрицы У2Ед. Так как Н = вV2Eд, то Л, = в^з, а следовательно,

Отсюда

ЛЛ, Л,

"Л! = = в

Лв >п|А = Лв £ ^ + а) = 1 £ ЛТ^.

,=1 ,=1

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

2вЕ

МР

N —

Л,

^ Л, + а

, = 1 и

N -

(15)

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

3. Процедура поиска оптимальной модели

Поиск оптимальной модели происходит на множестве порождаемых моделей на каждой итерации алгоритма. Перед работой алгоритма заданы множество измеряемых данных О и множество гладких функций С. Задан начальный набор конкурирующих моделей ^о = {/1,..., /м|/ £ Ф}, в котором каждая модель /г есть суперпозиция функций {дгз-}Г==1. Каждой функции дгз- — элементу модели /г ставится в соответствие гиперпараметр агз-, характеризующий начальную плотность распределения вектора параметров Ъгз- этой функции. Каждой модели /г поставлен в соответствие гиперпараметр вг начального приближения. Параметры начального приближения для г-й модели назначаются исходя из априорного распределения данных, определяемых значением вг. Далее выполняется последовательность шагов, приведенных ниже, которые повторяются заданное количество раз.

1. Методом сопряженных градиентов [9] минимизируются штрафные функции для каждой модели /г, г = 1,..., М. Отыскиваются параметры моделей wMP.

2. После нахождения параметров wMP исходя из (14) и (15) определяются новые значения гиперпараметров — аПзеш и вгпеш. Гиперпараметр вг функции /г вычисляется для всего набора данных и равен

N - Ъ

пеш

в (/г)"

Гиперпараметр агз- вычисляется для каждой функции дгз- из суперпозиции /г и равен

апеш =

™ - 1г

%3 (Ъ,- ) .

Здесь значения функционалов 7г и (Ъгз-) вычисляются только для подмножества тех параметров Ъгз- из множества wi, которые являются параметрами функции дгз-. Изменение гиперпараметров повторяется итерационно до тех пор, пока локальный минимум не

останется постоянным.

3. Заданы следующие правила построения производных моделей /1,..., /М. Для каждой модели /г строится производная модель /г'. В /г выбирается функция дгз- с наименьшим значением агз-. Выбираются произвольная модель / из Е0 \{/г} и ее произвольная функция д^. Модель /' порождается из модели / путем замещения функции дгз- с ее аргументами на функцию д^ с ее аргументами.

4. С заданной вероятностью п каждая модель /' подвергается изменениям. В изменяемой модели выбирается ^'-я функция, причем закон распределения вероятности выбора функции ) задан. Из множества С случайным образом выбирается функция д', которая замещает функцию дз-. Гиперпараметр агз- этой функции определяется как шах(агз-).

з

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

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

4. Численный эксперимент

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

0.9

measure data regression

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

-0.8 -0.6 -0.4 -0.2

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

показаны исходные данные, штриховой — значения модели 2. По оси абсцисс отложено значение свободной переменной, по оси ординат — значение зависимой переменной. Выборка, представленная данной кривой, содержит четыре тысячи отсчетов. Для верификации полученных моделей использовалось 118 выборок.

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

Таблица 1. Множество С базовых функций

№ Функция Описание Параметр

Функции двух переменных аргументов, д(Ь,Х1,Х2)

1 plus y = xi + x2 —

2 times y = xix2 —

3 divide y = xi/x2 —

Функции одного переменного аргумента, д(Ъ,х\)

7

8

9

10

multiply add

gaussian

linear

parabolic

cubic

logsig

(x - o2

2a2

y = ax y = x + a Л

y = - exP

V 2na

y = ax + b

y = ax2 + bx + c

y = ax3 + bx2 + cx + d

=_л_

y = 1+exp(-a(x - 0)

+ a

+ a

Л,a, a

a, b a, b, c a, b, c, d

Л,a, a

6

Таблица 2. Описание выбранных моделей

Описание Модель

1 2 3

Ошибка р1 0.0034 0.0037 0.0035

Ошибка р2 0.0421 0.0325 0.00338

Число параметров 16 16 16

Структура модели + х ±

Н + х + Н Н Н х + X + х 1 ± Н х Н Н х х Н + х X X Н Нх х х х

Примечание: Н — gaussian, с — еиЫе, I — Нпеаг, +--р1ш,

х — times, X--divide.

Выбор моделей производился из более тысячи порожденных моделей. В табл. 2 приведены три модели, полученные в результате работы алгоритма. Качество моделей оценивалось по ошибкам р1,р2 и числу параметров в векторе параметров Значения ошибок каждой модели получены путем усреднения результатов оптимально настроенной модели по 118 выборкам. Ошибка р1 — среднеквадратичная относительная ошибка

Р1

\

(у - / (х.)

N ^ V тах(У.)

ошибка р2 — максимальная относительная ошибка

- /(х.)|

р2 = тах

.=1.....N тах(у.)

В качестве примера рассмотрим модель 2. Она состоит из суперпозиции восьми функций /2 = д1(д2(д3(д4(д5(х),д6(х)),д7(х)),х),д8(х)). Функции сложения и умножения д1 = х(0, •, •) и д2, ...,д4 = + (0, •, •) имеют первым аргументом пустой вектор параметров;

а. ( (х - е.)2~

д5,...,д7 = Н(Ь., •), г = 1,..., 3, и д8 = /(Ь4, •). Функции Н = ,—— ехр I--—2— ) имеют

\j2nui \ 2а, ) векторы параметров Ь. = (А.,ц.,а.,а.), а функция I = (ах + Ь) имеет вектор параметров Ь4 = (а, Ь).

Модель /2 можно переписать в виде

/(', х) = /(Ь4,х) 1 х х + ^ Н(Ь.,х)

V .=1

где х = х и ' = Ь1.Ь2.Ь3.Ь4. Развернутый вид модели

< , и\-1 ( , ^ А. / (х - С.)2\ , у = (ах+Ь) ^+£ ^ех^+а

Модель /2 была использована экспертами для анализа и прогноза концентрации кислорода в выхлопных газах дизельного двигателя.

2

3

Заключение

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

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

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

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

Список литературы

[1] MALADA H.R., Ivakhnenko A.G. Inductive Learning Algorithms for Complex Systems Modeling. CRC Press, 1994.

[2] Bishop C.M., Tipping M.E. Bayesian regression and classification // Advances in Learning Theory: Methods, Models and Applications / J. Suykens, G. Horvath et. al. (Eds). IOS Press, NATO Sci. Ser. III: Computer and Systems Sciences, 2000. Vol. 190. P. 267-285.

[3] LeCun Y., Denker J.S., Solla S.A. Optimal brain damage // Advances in Neural Information Processing Systems / D.S. Touretzky (Ed.). Morgan Kaufmann, San Mateo, CA, 1990. P. 598-605.

[4] MaoKay D. Information, Inference, Learning Algorithms. Cambridge: Cambridge Univ. Press, 2003.

[5] MaoKay D. Hyperparameters: optimise or integrate out? // Maximum Entropy and Bayesian Methods / G. Heidberger (Ed.). Santa Barbara; Dordrecht: Kluwer, 1993.

[6] MaoKay D. Bayesian interpolation // Neural Comp. 1992. Vol. 4, N 3. P. 415-447.

[7] Nabney I.T. NETLAB: Algorithms for Pattern Recognition. N.Y.; Berlin: Springer-Verl., 2004. P. 330.

[8] MaoKay D. Choice of basis for Laplace approximation // Machine Learning. 1998. Vol. 33, N 1.

[9] Branch M.A., Coleman T.F., Li Y. A subspace, interior, and conjugate gradient method for large-scale bound-constrained minimization problems // SIAM J. on Sci. Comp. 1999. Vol. 21, N 1. P. 1-23.

Поступила в редакцию 27 июля 2006 г.