\year{2005} \title{Поиск оптимальной регрессионной модели в индуктивно заданном множестве}% \thanks{Работа выполнена при финансовой поддержке РФФИ (грант \mbox{№\,04-01-00401}).} \authors{В.~В.~Стрижов, канд.~физ.-мат.~наук\\ (Вычислительный центр им. А.~А. Дородницына РАН)} \abstract{% Описана процедура поиска оптимальной параметрической регрессионной модели в классе моделей, определенном суперпозициями гладких функций из заданного множества. Для поиска используются оценки плотности распределения параметров элементов моделей. Параметры моделей оцениваются с помощью методов нелинейной оптимизации. Для иллюстрации приведена задача о моделировании изменения давления в камере внутреннего сгорания дизельного двигателя.} \begin{document} \maketitle \section{Введение} Проблема выбора оптимальной регрессионной модели имеет большую историю, однако остается одной их самых актуальных в области распознавания образов. А.Г.~Ивахненко, еще в 1968 году, предложил метод группового учета аргументов, МГУА~\cite{ivakhnenko94}. Согласно этому методу модель, доставляющая наилучшее приближение, отыскивается во множестве последовательно порождаемых моделей. В~частности, для построения моделей как суперпозиций функций, использовались полиномиальные функции, нейронные сети и некоторые другие функции. А.Г.~Ивахненко и его ученики создали ряд алгоритмов синтеза моделей и предложили методы оценки качества моделей. При синтезе конкурирующих моделей появляется задача определения значимости элементов модели. В работе К.~Бишопа~\cite{bishop03} предложен метод анализа распределения параметров однослойных нейронных сетей посредством гиперпараметров, то есть параметров распределения параметров аппроксимирующих функций. Для каждого элемента нейронной сети оценивается плотность Гауссовского распределения его параметров и делается вывод о том, насколько информативен данный элемент исследуемой регрессионной модели. Для модификации моделей Ле Кюн предложил метод, называемый методом оптимального отсечения (optimal brain damage)~\cite{lecun90}. Этот метод состоит в исключении некоторых, наименее информативных, элементов регрессионной модели с тем условием, что при таком исключении качество аппроксимации уменьшается незначительно. При исключении отдельных элементов модели становится возможным оценить вклад этих элементов по значениям заданной функции качества аппроксимации. %Проблема выбора свободных переменных на каждом элементе~$g$ %суперпозиции является общей проблемой для задач построения %регрессионных моделей. (((В контексте данной работы свободными %переменными также являются значения элементов на каждом уровне %вложенности суперпозиции.))) Для линейных моделей существует %хорошо разработанная процедура анализа главных компонент, в %которой наиболее важные переменные имеют б\'ольшую корреляцию с %первой главной компонентой. В нелинейных моделях такой подход не %может быть использован. Те свободные переменные, которые не имеют %никакого влияния на зависимую переменную, могут иметь произвольные %веса. %Более удачный %метод, предложенный Ле~Кюн~\cite{lecun90} для нейронных сетей, %заключается в анализе матрицы Гессе и использует алгоритм %сокращения весов. При Байесовском подходе мы можем связать %отдельные гиперпараметры с группами весов. Каждый гиперпараметр %будет представлять единицу деленную на дисперсию весов данного %входа. В течение Байесовского обучения есть возможность изменять %гиперпараметры. Так как гиперпараметр есть обратная величина %дисперсии, по его малому значению мы можем понять, что большие %веса допустимы и сделать вывод о важности данной свободной %переменной. Это оставляет открытым вопрос, какие веса можно %считать больш\'ими, но уже можно сказать, какие переменные более %важны. Проблема сравнения и выбора регрессионных моделей получила новое развитие после ряда публикаций Д.~МакКая~[4---6] %\cite{mackay,mackay94,mackay98,mcakay92a} предложившего при выборе модели из заданного множества использовать не информационные критерии, например AIC~--- Akakie Information Criterion, а двухуровневый Байесовский вывод и правило Оккама. На первом уровне вывода вычисляются плотности вероятностей распределения параметров каждой модели из заданного множества. На втором уровне вывода вычисляется правдоподобие моделей. Правило Оккама состоит в том, что вероятность выбора более сложной модели меньше, чем вероятность выбора более простой модели при сравнимом значении функции качества аппроксимации. Предлагаемый метод заключается в следующем. Поиск моделей выполняется по итерационной схеме ``порождение-выбор'' в соответствии с определенными правилами порождения моделей и критерием выбора моделей. Последовательно порождаются наборы конкурирующих моделей. Каждая модель в наборе является суперпозицией элементов заданного множества гладких параметрических функций. После построения модели, каждому элементу суперпозиции ставится в соответствие гиперпараметр. Параметры и гиперпараметры модели последовательно настраиваются. Из набора выбираются наилучшие модели для последующей модификации. При модификации моделей, по значениям гиперпараметров делаются выводы о целесообразности включения того или иного элемента в модель следующего порождаемого набора. %Алгоритм, описанный ниже, построен на основе метода %группового учета аргументов Ивахненко~\cite{ivakhnenko94} и %Байесовского подхода в нейронных сетях с анализом весовых %коэффициентов, описанного в %статье Бишопа~\cite{bishop03}. %Наиболее полная подборка работ по %данной тематике опубликована на сайтах~\cite{wwwGMDH, wwwMacKay} %Наиболее близкой публикацией является~\cite{nikolaev}. Поставим задачу нахождения регрессионной модели нескольких свободных переменных следующим образом. Задана выборка~--- множество $\{\x_1,...,\x_N|\x\in\R^M\}$ значений свободных переменных и множество $\{y_1,...,y_N| y\in\R\}$ соответствующих им значений зависимой переменной. Обозначим оба эти множества как множество исходных данных~$D$. Также задано множество $G=\{g|g:\R\times...\times\R\longrightarrow\R\}$ параметрических гладких функций $g=g(\be,\cdot,\cdot,...,\cdot)$. Первый аргумент функции~$g$~--- вектор-строка параметров~$\be$, последующие~--- переменные из множества действительных чисел, рассматриваемые как элементы вектора свободных переменных. Рассмотрим произвольную суперпозицию, состоящую из не более чем~$r$ функций~$g$. Эта суперпозиция задает параметрическую регрессионную модель~% ~$f$ %\begin{equation}\label{eq0} $f=f(\w,\x).$ %\end{equation} Модель~$f$ зависит от вектора свободных переменных~$\x$ и от вектора параметров~$\w$. Вектор~$\w\in\R^W$ состоит из присоединенных векторов-параметров функций~$g_1,...,g_r$, то есть, $\w=\be_1\vdots\be_2\vdots...\vdots\be_r$, где $\vdots$~--- знак присоединения векторов. Обозначим~$\Phi=\{f_i\}$~--- множество всех суперпозиций, индуктивно порожденное элементами множества~$G$. Требуется найти такую модель~$f_i$, которая доставляет максимум функционала~$p(\w|D,\alpha,\beta,f_i)$. Этот функционал, определяемый далее, включает искомую модель~$f_i(\w,\x)$ и ее дополнительные параметры~$\alpha$ и~$\beta$. \section{Выбор регрессионных моделей и гипотеза порождения данных} Общий подход к сравнению нелинейных моделей описан МакКаем~\cite{mackay} и заключается в следующем. Рассмотрим набор конкурирующих моделей $f_1,...,f_M$. Априорная вероятность модели~$f_i$ определена как $P(f_i)$. При появлении данных~$D$ апостериорная вероятность модели~$P(f_i|D)$ может быть найдена по теореме Байеса, %$$P(f_i|D)=\frac{P(f_i)p(D|f_i)}{p(D)},$$ $$P(f_i|D)=\frac{P(f_i)p(D|f_i)}{\sum_{i=1}^Mp(D|f_i)P(f_i)},$$ где~$p(D|f_i)$~--- функция соответствия модели данным. Знаменатель дроби обеспечивает выполнение условия~$\sum_{i=1}^MP(f_i|D)=1$. Вероятности моделей~$f_1$ и~$f_2$, параметры которых идентифицированы по данным~$D$, сравнимы как \begin{equation} \frac{P(f_1|D)}{P(f_2|D)}=\frac{P(f_1)p(D|f_1)}{P(f_2)p(D|f_2)}. \label{bayes2}% \end{equation}% В данном отношении не указаны вероятности появления данных~$D$, поскольку для обеих моделей эти вероятности равны. Отношение~$\frac{p(D|f_1)}{p(D|f_2)}$~--- есть отношение правдоподобия моделей. Отношение~$\frac{P(f_1)}{P(f_2)}$ является априорной оценкой предпочтения одной модели другой. При моделировании отдается предпочтение наиболее простым и устойчивым моделям. Но если априорные оценки~$P(f_i)$ моделей одинаковы, то есть, нет причины предпочитать одну модель другой, %(см.~\cite{mcakay92a}) то их необходимо сравнивать по значениям $p(D|f_i)$: %\begin{equation} $$ p(D|f_i)=\int{}p(D|\w,f_i)p(\w|f_i)d\w. %\label{intpd} $$%\end{equation} % Апостериорная плотность распределения параметров~$\w$ функции~$f_i$ при заданной выборке~$D$ равна \begin{equation} p(\w|D,f_i)=\frac{p(D|\w,f_i)p(\w|f_i)}{p(D|f_i)}, \label{bayes1}% \end{equation}% где~$p(\w|f_i)$~--- априорно заданная плотность вероятности параметров начального приближения, $p(D|\w,f_i)$~--- функция правдоподобия параметров модели, а знаменатель~$p(D|f_i)$ обеспечивает выполнение условия~$\int{}p(\w|D,f_i)d\w=1$. Он задан интегралом в пространстве параметров $\int{}p(\w'|D,f_i)p(\w'|f_i)d\w'.$ Формулы~(\ref{bayes1}) и~(\ref{bayes2}) называются формулами Байесовского вывода первого и~второго уровня. Рассмотрим регрессию $y=f_i(\be,\x)+\nu$ с аддитивным Гауссовским шумом с дисперсией~$\sigma_\nu$ и с нулевым матожиданием. Тогда плотность вероятности появления данных % %\begin{equation} $$ p(y|x,\w,\beta,f_i)\triangleq{}p(D|\w,\beta,f)=\frac{\exp(-\beta{}E_D(D|\w,f_i))}{Z_D(\beta)}, $$ %\label{pdw} %\end{equation} где $\beta=\frac{1}{\sigma^2_\nu}$. Нормирующий множитель~$Z_D(\beta)$ задан выражением \begin{equation} Z_D(\beta)=\left(\frac{2\pi}{\beta}\right)^{\frac{N}{2}} \label{ZD} \end{equation} и взвешенный функционал ошибки в пространстве данных \begin{equation} \beta{}E_D=\frac{\beta}{2}\sum_{n=1}^N(f_i(\x_n)-y_n)^2. \label{bED} \end{equation} %Выражение~(\ref{pdw}) в развернутой записи имеет вид %$$ %p(D|\w,\beta,f)=\frac{1}{(\sigma_\nu\sqrt{2\pi})^N}\exp\left(-\frac{\sum\limits_{n=1}^N(f(\x_n)-y_n)^2}{2\sigma_\nu^2}\right). %$$ %Задача оценки наиболее правдоподобных весов в сложных моделях не %является корректной в смысле Адамара. Введем регуляризующий параметр~$\alpha$, который отвечает за то, насколько хорошо модель должна соответствовать зашумленным данным. Функция плотности вероятности параметров с заданным гиперпараметром~$\alpha$ имеет вид $$ p(\w|\alpha,f_i)=\frac{\exp(-\alpha{}E_W(\w|f_i))}{Z_W(\alpha)}, $$ где~$\alpha$~--- обратная дисперсии распределения параметров, $\alpha=\sigma_\w^{-2}$, а нормирующая константа~$Z_W$ определена дисперсией распределения параметров как \begin{equation} Z_W(\alpha)=\left(\frac{2\pi}{\alpha}\right)^\frac{W}{2}. \label{ZW} \end{equation} Требование к малым значениям параметров~\cite{netlab} предполагает Гауссовское априорное распределение с нулевым средним: $$ p(\w)=\frac{1}{Z_W}\exp{(-\frac{\alpha}{2}||\w||^2)}. $$ Так как переменные~$\alpha$ и~$\beta$ являются параметрами распределения параметров модели, в дальнейшем будем называть их гиперпараметрами. Исключая нормирующую константу~$Z_W$, которая не зависит от параметров~$\w$ и логарифмируя, получаем \begin{equation} \alpha{E}_W=\frac{\alpha}{2}||\w||^2.% \label{aED} \end{equation} Эта ошибка регуляризирует параметры, начисляя штраф за их чрезмерно большие значения. %Функционал $E_W$ называется %регуляризирующим~\cite{gull}. При заданных значениях гиперпараметров~$\alpha$ и~$\beta$ выражение~(\ref{bayes1}) для фиксированной функции~$f_i$ будет иметь вид $$ p(\w|D,\alpha,\beta)=\frac{p(D|\w,\beta)p(\w|\alpha)}{p(D|\alpha,\beta)}. $$ Записывая функцию ошибки в виде $S(\w)=\alpha{}E_W+\beta{}E_D$, получаем \begin{equation} p(\w|D,\alpha,\beta,f_i)=\frac{\exp(-S(\w|f_i))}{Z_S(\alpha,\beta)}, \label{pwabd} \end{equation} где $Z_S$~--- нормирующий множитель. % \section{Нахождение параметров модели} Рассмотрим итеративный алгоритм для определения оптимальных параметров~$\w$ и гиперпараметров~$\alpha,\beta$ при заданной модели~$f_i$. %При этом интегрирования по всему пространству параметров выполняться %не будет и, следовательно, этот алгоритм не будет полноценными %Байесовским подходом. % see type II maximum likelyhood method, Berger 1985, p.341 netlab Корректный подход заключается в интегрировании всех неизвестных параметров и гиперпараметров. Апостериорное распределение параметров определяется как \begin{equation} p(\w|D)=\iint{}p(\w,\alpha,\beta|D)d\alpha{}d\beta=\iint{}p(\w|\alpha,\beta,D)p(\alpha,\beta|D)d\alpha{}d\beta, \label{intalpbet} \end{equation} что требует выполнить интегрирование апостериорного распределения параметров~$p(\w|\alpha,\beta,D)$ по пространству, размерность которого равна количеству параметров. Вычислительная сложность этого интегрирования весьма велика. Интеграл может быть упрощен при подходящем выборе начальных % точнее: априорных значений гиперпараметров: $p(\alpha)=1/\alpha$. %Тогда %интеграл может быть вычислен аналитически~\cite{mackay94}. %\cite{wolpert93} Приближение интеграла заключается в том, что апостериорная плотность распределения гиперпараметров~$p(\alpha,\beta|D)$ имеет выраженный пик в окрестности наиболее правдоподобных значений гиперпараметров~$\alpha^\m,\beta^\m$. Это приближение известно как аппроксимация Лапласа~\cite{mackay98}. При таком допущении интеграл~(\ref{intalpbet}) упрощается до $$ p(\w|D)\approx{}p(\w|\alpha^\m,\beta^\m,D)\iint{p(\alpha,\beta|D)}d\alpha{}d\beta\approx{}p(\w|\alpha^\m,\beta^\m,D). $$ Необходимо найти значения гиперпараметров, которые оптимизируют апостериорную плотность вероятности параметров, а затем выполнить все остальные расчеты, включающие~$p(\w|D)$ при фиксированных значениях гиперпараметров. % У МакКая была работа о том, почему их нельзя оптимизировать совместно Для нахождения функционала~$p(\w|\alpha,\beta,D)$, который использует апостериорное распределение параметров, рассмотрим аппроксимацию ошибки~$S(\w)$ на основе рядов Тейлора второго порядка: \begin{equation} S(\w)\approx{}S(\w^\m)+\frac{1}{2}(\w-\w^\m)^TA(\w-\w^\m). \label{approxs} \end{equation} В выражении~(\ref{approxs}) нет слагаемого первого порядка, так как предполагается, что~$\w^\m$ определяет локальный минимум функции ошибки, то есть $$\frac{\partial{}S(\w^\m)}{\partial{}w_\xi}=0$$ для всех значений~$\xi$. Матрица~$A$~--- это матрица Гессе функции ошибок: $$ A=\nabla^2{}S(\w^\m)=\beta\nabla^2{}E_D(\w^\m)+\alpha{}I. $$ %где элемент~$a_{ij}$ матрицы~$A$ вторых производных имеет вид %$$a_{ij}=\frac{\partial^2 S}{\partial{}w_i\partial{}w_j}.$$ Обозначим первое слагаемое правой части через~$H$. Подставив полученное приближенное значение $S(\w)$ в~(\ref{pwabd}) и обозначив~$\Delta\w=\w-\w^\m$, получим $$ p(\w|\alpha,\beta,D)=\frac{1}{\hat{Z}_S}\exp\left(-S(\w^\m)-\frac{1}{2}\Delta\w^TA\Delta\w\right), $$ Оценим нормирующую константу~$\hat{Z}_S$, необходимую для аппроксимации кривой Гаусса, как \begin{equation} \hat{Z}_S=\exp(-S(\w^\m))(2\pi)^\frac{W}{2}(\det{}A)^{-\frac{1}{2}}. \label{ZS} \end{equation} %Так как~(\ref{ZS}) не зависит от искомых пиков параметров~$\alpha$ %и~$\beta$, мы можем его опустить. Максимизируем функцию~$p(D|\alpha,\beta)$, изменяя значения гиперпараметров~$\alpha$ и~$\beta$. Это можно выполнить, интегрируя функцию плотности вероятности данных по пространству параметров~$\w$: \begin{equation} p(D|\alpha,\beta)=\int{}p(D|\w,\alpha,\beta)p(\w|\alpha,\beta)d\w=\int{}p(D|\w,\alpha,\beta)p(\w|\alpha)d\w, \label{intDab} \end{equation} где второй интеграл справедлив по причине того, что распределение параметров не зависит от дисперсии шума в силу гипотезы о Гауссовском распределении шума. %остается проблема назначения априорной плотности вероятности~$p(\alpha,\beta)$. Для упрощения вычислений мы допускаем, что распределение~$p(\alpha,\beta)$ является равномерным. %и поэтому может быть исключена из последующего анализа. Используя~(\ref{bED}),~(\ref{aED}), запишем~(\ref{intDab}) в виде $$ p(D|\alpha,\beta)=\frac{1}{Z_D(\beta)}\frac{1}{Z_D(\alpha)}\int\exp(-S(\w)))d\w. $$ Из ~(\ref{ZD}),~(\ref{ZW}),~(\ref{ZS}) и предыдущего выражения получим \begin{equation} \ln{}p(D|\w,\alpha,\beta)=-\alpha{}E_W^\m-\beta{}E_D^\m-\frac{1}{2}\ln|A|+\frac{W}{2}\ln{\alpha}+\frac{N}{2}\ln\beta-\frac{N}{2}\ln{}(2\pi). \label{lnpDwab} \end{equation} Для того, чтобы оптимизировать это выражение относительно~$\alpha$, найдем производную \begin{equation} \frac{d}{d\alpha}\ln|A|=\frac{d}{d\alpha}\ln\left(\prod_{j=1}^W\lambda_j+\alpha\right)% =\frac{d}{d\alpha}\sum_{j=1}^W\ln(\lambda_j+\alpha)=\sum_{j=1}^W\frac{1}{\lambda_j+\alpha}=\mbox{tr}(A^{-1}). \label{dda} \end{equation} В этом выражении~$\lambda_1,...,\lambda_W$~--- собственные значения матрицы~$H$. Так как функция ошибки на данных не является квадратичной функцией параметров, как при линейной или RBF регрессии, то непосредственно оптимизировать величину $\alpha$ невозможно, гессиан~$Н$ не является константой, а зависит от параметров~$\w$. Так как мы принимаем $A=H+\alpha{}I$ для вектора~$\w^\m$, который зависит от выбора~$\alpha$, то собственные значения~$H$ косвенным образом зависят от~$\alpha$. Таким образом, формула~(\ref{dda}) игнорирует переменные в выражении $d\lambda_j/\alpha$. % С использованием этого приближения, производная~(\ref{lnpDwab}) с учетом~$\alpha$ равна $$ \ln{}p(D|\w,\alpha,\beta)=-E_W^\m-\frac{1}{2}\sum_{j=1}^W\frac{1}{\lambda_j+\alpha}+\frac{W}{2\alpha}. $$ Приравнивая последнее выражение к нулю и преобразовывая, получаем выражение для~$\alpha$ \begin{equation} 2\alpha{}E_W^\m=W-\sum_{j=1}^W\frac{\alpha}{\lambda_j+\alpha}. \label{alpha} \end{equation} Обозначим вычитаемое правой части через~$\gamma$ $$ \gamma=\sum_{j=1}^W\frac{\alpha}{\lambda_j+\alpha}. $$ Те компоненты суммы, в которых~$\lambda_j\gg\alpha$ привносят вклад, близкий к единице, а те компоненты суммы, в которых $0<\lambda_j\ll\alpha$, привносят вклад, близкий к нулю. %Следовательно, переменная~$\gamma$ является мерой числа хорошо %определенных параметров. % см. Бишоп, 1995, section 4. Для нахождения гиперпараметра~$\beta$ рассмотрим задачу оптимизации~(\ref{lnpDwab}). Обозначим через~$\mu_j$ собственные значение матрицы $\nabla^2{}E_D$. Так как $H=\beta\nabla^2{}E_D$, то $\lambda_j=\beta\mu_j$ и следовательно, $$ \frac{d\lambda_j}{d\beta}=\mu_j=\frac{\lambda_j}{\beta}. $$ Отсюда, $$ \frac{d}{d\beta}\ln|A|=\frac{d}{d\beta}\sum_{j=1}^W\ln(\lambda_j+\alpha)=\frac{1}{\beta}\sum_{j=1}^W\frac{\lambda_j}{\lambda_j+\alpha}. $$ Дифференцируя, как и в случае нахождения~$\alpha$, мы находим, что оптимальное значение~$\beta$ определено как \begin{equation} 2\beta{}E_D^\m=N-\sum_{j=1}^W\frac{\lambda_j}{\lambda_j+\alpha}=N-\gamma. \label{beta} \end{equation} Способ вычисления оптимальных значений гиперпараметров~$\alpha$ и~$\beta$ описан в следующем разделе. % \section{Процедура поиска оптимальной модели} Поиск оптимальной модели происходит на множестве порождаемых моделей на каждой итерации алгоритма. Перед работой алгоритма заданы множество измеряемых данных~$D$ и множество гладких % на области определения свободной переменой функций~$G$. Задан начальный набор конкурирующих моделей, $F_0=\{f_1,...,f_{M}|f\in\Phi\}$, в котором каждая модель~$f_i$ есть суперпозиция функций~$\{g_{ij}\}_{j=1}^{r_i}$. Каждой функции~$g_{ij}$~--- элементу модели~$f_i$ ставится в соответствие гиперпараметр~$\alpha_{ij}$, характеризующий начальную плотность распределения вектора параметров~$\be_{ij}$ этой функции. Каждой модели~$f_i$ поставлен в соответствие гиперпараметр~$\beta_{i}$ начального приближения. Параметры~$i$-й модели назначаются исходя из априорного распределения данных, определяемых значением~$\beta_i$. Далее выполняется последовательность шагов, приведенных ниже, которые повторяются заданное количество раз. 1. Методом сопряженных градиентов минимизируются штрафные функции $S_i(\w)$ для каждой модели~$f_i, i=1,...,M$~\cite{branch99}. Отыскиваются параметры моделей~$\w^\m_i$. % для заданных значений $\alpha_{i\cdot},\beta_i$. 2. После нахождения параметров~$\w^\m_i$ определяются, исходя из~(\ref{alpha}) и~(\ref{beta}), новые значения гиперпараметров~$\alpha_{ij}^\new$ и~$\beta_{i}^\new$. Гиперпараметр~$\beta_i$ функции~$f_i$ вычисляется для всего набора данных и равен % % %$\beta_i=\frac{1}{2}(N-\gamma)E_{D}^{-1}(f_i).$ $$\beta_i^\new=\frac{N-\gamma_i}{E_{D}(f_i)}.$$ % Гиперпараметр~$\alpha_{ij}$ вычисляется для каждой функции~$g_{ij}$ из суперпозиции~$f_i$ и равен % % %$\alpha_{ij}=\frac{1}{2}\gamma_{ij}{}E_W(\be_{ij})^{-1}.$ $$\alpha_{ij}^\new=\frac{\gamma_i}{E_W(\be_{ij})}.$$ % Здесь значения функционалов~$\gamma_i$ и~$E_W(\be_{ij})$ вычисляется только для подмножества тех параметров~$\be_{ij}$ из множества~$\w_i$, которые являются параметрами функции~$g_{ij}$. Изменение гиперпараметров повторяется итерационно до тех пор пока локальный минимум~$S_i(\w)$ не останется постоянным. %После отыскания параметров и гиперпараметров моделей выполняются %действия генетического оптимизационного алгоритма: кроссовер, %мутация и селекция. Модификация классического %алгоритма~\cite{Goldberg89} заключается в следующем. 3. Заданы следующие правила построения производных моделей $f'_1,...,f'_{M}$. Для каждой модели~$f_i$ строится производная модель~$f'_i$. В~$f_i$ выбирается функция~$g_{ij}$ с наименьшим значением~$\alpha_{ij}$. Выбирается произвольная модель~$f_\xi$ из $F_0\backslash\{f_i\}$ и ее произвольная функция~$g_{\xi\zeta}$. Модель~$f'$ порождается %производится из модели~$f$ путем замещения функции~$g_{ij}$ с ее аргументами на функцию~$g_{\xi\zeta}$ с ее аргументами. \begin{figure}[!bp] \centering{\includegraphics[angle=0,width=0.6\textwidth]{./model145.eps}} \caption{Исходная выборка и восстановленная выборка, полученная моделью \No~2}% \label{fig1} \end{figure} 4. С заданной вероятностью~$\eta$ каждая модель~$f'_i$ подвергается изменениям. В изменяемой модели выбирается~$j$-тая функция, причем закон распределения вероятности выбора функции~$p(j)$ задан. Из множества~$G$ случайным образом выбирается функция~$g'$ и замещает функцию~$g_j$. Гиперпараметр~$\alpha_{ij}$ этой функции определяется как $\max\limits_j(\alpha_{ij})$. Вектор параметров этой функции~$\be_{ij}$ равен нулю или назначается при задании~$G$. 5. При выборе моделей из объединенного множества родительских и порожденных моделей в соответствии с критерием~$S(\w)$ выбираются~$M$ наилучших, которые используются в дальнейших итерациях. % % \section{Численный эксперимент} Ниже описывается пример построения регрессионной модели. Объектом моделирования является кривая одной свободной переменой, представленная набором измерений давления в камере внутреннего сгорания дизельного двигателя. На рис.~\ref{fig1} сплошной кривой показаны исходные данные, пунктиром показаны значения модели~\No~2. По оси абсцисс отложено значение свободной переменной, по оси ординат~--- значение зависимой переменной. Выборка, представленная данной кривой, содержит четыре тысячи отсчетов. Для верификации полученных моделей использовалось 118 выборок. Экспертами было задано множество базовых функций~$G$ из элементов которого порождаются регрессионные модели. Список функций приведен в таблице~\ref{tbl1}. Множество~$F_0$ моделей начального приближения также было задано экспертами. % \begin{table}[!htbp] \centering{ % {\footnotesize \begin{tabular}{|c|l|l|l|} \hline % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ... \No & Функция & Описание & Параметры \\ \hline\hline \multicolumn{4}{|c|}{Функции двух переменных аргументов, $g(\be,x_1,x_2)$}\\ \hline 1 & plus & $y=x_1+x_2$ & --\\ %2 & minus & $y=x_1-x_2$ & --\\ 2 & times & $y=x_1x_2$ & --\\ 3 & divide & $y=x_1/x_2$ & --\\ \hline \multicolumn{4}{|c|}{Функции одного переменного аргумента, $g(\be,x_1)$}\\ \hline 4 & multiply & $y=ax$ & $a$\\ 5 & add & $y=x+a$ & $a$\\ 6 & gaussian & $y=\frac{\lambda}{\sqrt{2\pi}\sigma}\mbox{exp}\left({-\frac{(x-\xi)^2}{2\sigma^2}}\right)+a$ & $\lambda,\sigma,\xi,a$ \\% 7 & linear & $y=ax+b$ & $a,b$\\ 8 & parabolic & $y=ax^2+bx+c$ & $a,b,c$\\ 9 & cubic & $y=ax^3+bx^2+cx+d$ & $a,b,c,d$\\ 10& logsig & $y=\frac{\lambda}{1 + \exp(-\sigma(x-\xi))}+a$ & $\lambda,\sigma,\xi,a$\\ %10 & tansig & $y=\frac{\lambda}{1+\exp(-2\sigma(x-\xi))}+a$ & $\lambda,\sigma,\xi,a$\\ \hline \end{tabular}}% end of centering } %end of footnotesize \caption{Множество~$G$ базовых функций}% \label{tbl1} \end{table} Выбор моделей производился из более тысячи порожденных моделей. В таблице~\ref{tbl2} приведены три модели, полученные в результате работы алгоритма. Качество моделей оценивалось по ошибкам~$\rho_1, \rho_2$ и числу параметров в векторе параметров~$\w$. Значения ошибок каждой модели получены путем усреднения результатов оптимально настроенной модели по~118 выборкам. Ошибка~$\rho_1$~--- среднеквадратичная относительная ошибка $$ \rho_1=\sqrt{\frac{1}{N}\sum_{i=1}^n\left(\frac{y_i-f(\x_i)}{\max(y_i)}\right)^2}, $$ ошибка ~$\rho_2$~--- максимальная относительная ошибка $$ \rho_2=\max_{i=1,...,N}\frac{|y_i-f(\x_i)|}{\max(y_i)}. $$ \begin{table}[!htbp] \centering {\footnotesize% \drawwith{\dottedline{1}}% \setlength{\GapDepth}{1mm}% \setlength{\GapWidth}{2mm} \begin{tabular}{|p{3cm}|c|c|c|} \hline% \No~модели & 1 & 2 & 3\\ \hline% Ошибка $\rho_1$ & 0.0034 & 0.0037 & 0.0035\\ \hline% Ошибка $\rho_2$ & 0.0421 & 0.0325 & 0.00338\\ \hline% Число параметров & 16 & 16 & 16\\ \hline% Описание & %1 \begin{bundle}{$+$}\chunk{\begin{bundle}{$h$}\chunk{$x$}\end{bundle}}\chunk{\begin{bundle}{$+$}\chunk{\begin{bundle}{$+$}\chunk{\begin{bundle}{$h$}\chunk{$x$}\end{bundle}}\chunk{\begin{bundle}{$h$}\chunk{$x$}\end{bundle}}\end{bundle}}\chunk{\begin{bundle}{$h$}\chunk{$x$}\end{bundle}}\end{bundle}}\end{bundle} % 2 & \begin{bundle}{$\times$}\chunk{\begin{bundle}{$+$}\chunk{\begin{bundle}{$+$}\chunk{\begin{bundle}{$+$}\chunk{\begin{bundle}{$h$}\chunk{$x$}\end{bundle}}\chunk{\begin{bundle}{$h$}\chunk{$x$}\end{bundle}}\end{bundle}}\chunk{\begin{bundle}{$h$}\chunk{$x$}\end{bundle}}\end{bundle}}\chunk{$x$}\end{bundle}}\chunk{\begin{bundle}{$\div$}\chunk{\begin{bundle}{$l$}\chunk{$x$}\end{bundle}}\end{bundle}}\end{bundle}% 3 & \begin{bundle}{$+$}\chunk{\begin{bundle}{$\times$}\chunk{\begin{bundle}{$+$}\chunk{\begin{bundle}{$\div$}\chunk{\begin{bundle}{$h$}\chunk{$x$}\end{bundle}}\chunk{\begin{bundle}{$c$}\chunk{$x$}\end{bundle}}\end{bundle}}\chunk{\begin{bundle}{$h$}\chunk{$x$}\end{bundle}}\end{bundle}}\chunk{$x$}\end{bundle}}\chunk{\begin{bundle}{$h$}\chunk{$x$}\end{bundle}}\end{bundle} \\ \hline % \multicolumn{4}{l} {\footnotesize Легедна: h~--- gaussian, c~--- cubic, l~--- logsig,} \\% \multicolumn{4}{l} {\footnotesize $+$~--- plus, $\times$~--- times,$\div$~--- divide} \\% % \end{tabular} }% \caption{Описание выбранных моделей} \label{tbl2} \end{table} В строке <<Описание>> таблицы~\ref{tbl2} показана структура модели в виде дерева. В качестве примера рассмотрим модель \No 2. Эта модель состоит из суперпозиции восьми функций~$f_2=g_1(g_2(g_3(g_4(g_5(x),g_6(x)),g_7(x)),x),g_8(x))$. Функции~$g_1=\times(\varnothing,\cdot,\cdot)$ и $g_2,...,g_4=+(\varnothing,\cdot,\cdot)$, сложения и умножения, имеют первым аргументом пустой вектор параметров; $g_5,...,g_7=h(\be_i,\cdot)$, $i=1,...,3$, и $g_8=l(\be_4,\cdot)$. Функции~$h=\frac{\lambda_i}{\sqrt{2\pi}\sigma_i}\mbox{exp}\left({-\frac{(x-\xi_i)^2}{2\sigma_i^2}}\right)$ имеют векторы параметров~$\be_i=\langle\lambda_i,\mu_i,\sigma_i,a_i\rangle$, и функция~$l=(ax+b)^{-1}$ имеет вектор параметров~$\be_4=\langle{a,b}\rangle$. Модель~$f_2$ можно переписать в виде $f(\w,\x)=l(\be_4,x)\times\left(x+\sum_{i=1}^3h(\be_i,x)\right),$ где~$\x=x$ и~$\w=\be_1\vdots\be_2\vdots\be_3\vdots\be_4$. % Развернутый вид модели: $$ y=(ax+b)^{-1}\left(x+\sum_{i=1}^3\frac{\lambda_i}{\sqrt{2\pi}\sigma_i}\mbox{exp}\left({-\frac{(x-\xi_i)^2}{2\sigma_i^2}}\right)+a_i\right). $$ Модель~$f_2$ была использована экспертами для анализа и прогноза концентрации кислорода в выхлопных газах дизельного двигателя. % \section*{Заключение} Универсальные регрессионные модели, например, нейронные сети или радиальные базисные функции, при обработке результатов измерений часто имеют большое число параметров и получаются пере\-обученными. Для достижения результатов в построении несложных и достаточно точных моделей была поставлена задача о выборе регрессионной модели, которая состоит из суперпозиции гладких функций. Для выбора наилучшей модели из индуктивно заданного множества использован двухуровневый Байесовский вывод. В связи со сложностью вычисления значений интегралов вывода, были предложены процедуры приближения, которые позволяют отыскивать адекватные модели за приемлемое время вычислений. Предложенная процедура выбора регрессионных моделей использует гиперпараметры, поставленные в соответствие элементам модели. Эти гиперпараметры указывают на важность элементов модели. На основе информации о важности элементов итерационно порождаются новые модели. Сложность моделей ограничивается автоматически при сравнении моделей. Описанный метод был протестирован на задаче по аппроксимации кривой, полученной в результате измерений давления в камере внутреннего сгорания дизельного двигателя. Получена модель с удовлетворительной погрешностью аппроксимации.