ИССЛЕДОВАНИЕ АЛГОРИТМОВ ВЫБОРА ОПТИМАЛЬНЫХ КООРДИНАТ УЗЛОВЫХ ТОЧЕК В ПОЛУПАРАМЕТРИЧЕСКИХ МОДЕЛЯХ ШТРАФНЫХ СПЛАЙНОВ Текст научной статьи по специальности «Математика»

Научный вестник НГТУ. - 2013. -№ 2(51) ОБРАБОТКА ИНФОРМАЦИИ

УДК 519.23

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

В.И. ДЕНИСОВ, B.C. ТИМОФЕЕВ, A.B. ФАДДЕЕНКОВ

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

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

ВВЕДЕНИЕ

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

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

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

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

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

Статья получена 27марта 2013 г.

Работа выполнена при финансовой поддержке РФФИ в рамках научного проекта №13-07-00299 а

1. сплаиновая полупараметрическая регрессия

Рассмотрим одну из известных полупараметрических регрессионных моделей следующего вида [3]:

yi = еЛ1+^+еА + р^ + р2^2 +...+р^* + , (1)

где у1 - значение отклика в г'-м наблюдении (г = 1,2,...,Д/ ); х^ - значение /-го регрессора в г-м наблюдении (у =1,2,...,т ); Р^...^ - неизвестные параметры;

+... + Втхы - параметрическая часть модели; Р1^1 + Р2^2 + ••• + Рк/гк ~ непараметрическая часть; /л , ■■■,//- значения базисных функций в г'-м наблюдении; г, - случайная ошибка в г'-м наблюдении (предполагается, что все ошибки независимы и имеют одинаковое

2 2

распределение с нулевым средним и дисперсией стЕ : б,- ~(0,сте), г -1,2,..., ТУ). В качестве базисных будем использовать функции

fw 2 /} (Х, ) 2 (Х, - bj У+ 2

| (Xj - bj )р, при (Xj ~bj)> О

npu(Xj-bj)< О

(2)

где bj е[а,й] -узловые точки (у = 1,2,...,к), р - некоторая положительная целая константа.

Непараметрическая часть модели (1) с базисными функциями (2) представляет собой сплайн порядка р, а саму модель в этом случае называют сплайновой регрессией.

В матричном виде уравнение (1) может быть представлено следующим образом:

У = X® + Е , (3)

где

X 2

41

*21

Х11

Y 2[ '2 4 ' Х1т (X1 - blY-i

' Х2т (Х2 - b1)4

Ум J .

Х1т (XN b1):+

(Х1 - bk )i (Х2 - bk

(Хм - bk

RR 2[e

0m ß1 4 ßk J >E2[rl r2

rN J •

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

в 2

=( X ux )-1

XX TY •

Естественно, что качество получаемой модели зависит от количества базисных точек и их координат. Однако рост числа узлов может приводить к излишней подгонке линии регрессии под исходные данные. Традиционным решением этой проблемы считается переход к так называемым «штрафным сплайнам» [3]. Идея этого метода заключается в том, что для снижения излишнего влияния непараметрической части, на ее параметры налагается ограничение (штраф) и вектор оценок параметров вычисляется следующим образом:

4

® 2(XTX + У2D

)-1 X UY,

(5)

где X - параметр сглаживания, О - (т + к)х(т + к) -матрица штрафа:

9 =

"0 • 0 0 • 0"

0 • 0 0 • 0

0 • 0 1

7 0

0 • 0 0 1.

0(тхт) 0(тхк)

'(кхк)

При X =0 сглаживание непараметрической части не проводится и оценка (5) совпадает

с обычной МНК-оценкой (4). Чрезмерное же увеличение параметра сглаживания (X —» +оо ) приводит к тому, что регрессионная модель (3) вырождается в модель, состоящую только из параметрической части. В связи с этим выбору величины параметра сглаживания следует уделять особое внимание.

Использование для этой цели метода максимального правдоподобия приводит к необходимости фиксации закона распределения случайных компонент модели (3) [7]. Наиболее часто предполагают нормальность распределения, однако на практике такое предположение далеко не всегда будет выполняться. В связи с этим в данной работе рассматриваются альтернативные способы оценивания, основанные на критериях кросс-валидации, обобщенной кросс-валидации, а также критерий Акаике.

Критерий кросс-валидации:

су = Х

1 =1

N Г ~ Л2 Ух - У

1 - Т

х,и у

(6)

где ц - г-й диагональный элемент матрицы = Х^ХТX + Х2И| ХТ. Критерий обобщенной кросс-валидации:

ССУ = -

N

Ъ(Уг - ')

г=1_

(1 - N(Ту ))

Критерий Акаике:

68 = 1о»{|(у,-у,)}*2^. (')

Определение оптимального значения для параметра сглаживания сводится к задаче одномерной минимизации статистик (6), (7) или (8). Отдельные результаты исследования качества восстановления регрессионной зависимости с использованием модели (3), (5) при различ-

2

ных X , к можно найти в [4].

Выбор координат узловых точек также оказывает существенное влияние на качество модели. Наиболее распространенным подходом для определения этих координат является следующий несложный алгоритм [3].

порядка

Алгоритм 1 - узловые точки bj ( ] — \,2,...,к ) определяются как выборочные квантили }

к + 1

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

деления наблюдений на отрезке [ Л ] и отсутствии повторных наблюдений координаты узловых точек могут быть определены следующим образом:

bw = а + (b - а)

к + 1

Очевидно, что в этом случае при выборе координат узловых точек не учитываются структура модели и имеющиеся значения зависимой переменной. Вопросы использования этой дополнительной информации применительно к модели (3-4) рассматривались в работе [5]. Предложенные ниже три алгоритма являются логическим обобщением подходов, рассмотренных в [5].

Алгоритм 2 - выбор узловых точек с использованием информационной матрицы модели (3). Начальное приближение для координат узловых точек bj ( j —1,2,.. ,,к) определяется

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

bw = Arg min det(M - 1 (j )),

bj£[a ,i]

где

м (j) 2 Ix u (j) X (j),

X (W) 2

'1 1 x2 1

XN1

'1 m (X1 - b1 >

('1 - bj )l

X2m (X2 b1)+

XNm (XN b1)

(X2 - bj )l

(XN - bj Y+

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

min det(M 1),

h ,---,bk е[ a>b]

(9)

1 _

где М 2 —X X - информационная матрица.

N

Алгоритм 3 - выбор узловых точек на основе остаточной суммы квадратов. Начальное приближение для координат / -й узловой точки (у —1,2,...,к) определяются как решение оптимизационной задачи

bf 2 Arg min ESS(j) ,

bw e[ a,b ]

где ESS (j) 2 e(j)Te(j), e(j) 2 Y - X(j) \ X(jf X(j)\' X(jfY.

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

min ESS, (10)

b\,4,bk е[ a,b]

где ESS 2 eUe, e = Y - X

x TX

X TY.

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

b w 2 Arg min det(D( j)),

bi£[a,b]

где D(j) 2 a£r (j) X(j)T X(j) , d^ (j) =

-1-1

ESS (j)

N - m - j

Дальнейшая оптимизация проводится как решение задачи

min det D,

b j",bk e[ h,b]

(11)

X rx

!

ESS

где 9 = дЕ , ,к

N - т - к

Комбинация различных методов выбора координат узловых точек и методов оценивания параметров модели (3) дает широкий спектр итоговых алгоритмов построения искомой линии регрессии.

2. ВЫЧИСЛИТЕЛЬНЫМ ЭКСПЕРИМЕНТ

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

Уг 2 У? +ri =Ро + ' - 1-8'г2 + +Pl( xi - ^1)+ + р2( 'v - h)+ • (12)

При моделировании отклика значения независимой переменной х равномерно варьировались на отрезке [0,1]. При каждой реализации набора исходных данных в функции (12) параметры , j — 0,1,2 определялись псевдослучайным датчиком как случайные величины, равномерно распределенные на отрезках [1,3], [0.25, 0.35], [0.6, 0.7] соответственно. Случайная

ошибка для каждого наблюдения генерировалась из предположений о нормальном распреде-

2 2 лении: '¿¡ ~ N(0, стЕ). Дисперсия ошибки стЕ выбиралась таким образом, чтобы величина

уровня шума

p=^L 100 %,

с

была равна наперед заданному значению. Уровень шума р введен в [ ] и определяется как отношение «шум»/«сигнал» в процентах, здесь стЕ - дисперсия ошибки s; 1 w

с —-(>'/ ~У ) \У i - не зашумленные измерения отклика).

На рис. 1 представлены примеры различных серий наблюдений, в каждой из которых параметры модели определялись по описанным выше правилам с уровнем шума р = 10 % .

Рис. 1. Примеры различных серий наблюдений, сгенерированных для модели (12)

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

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

ЕББ 2 (У - ¥)т (У - У), (13)

где У - вектор исходных данных, У - вектор оценок наблюдений, построенных по модели.

В качестве критерия точности соответствия истинной модели использовалась аналогичная сумма квадратов:

£55ист = (7ист-7)г(^ст-7), (14)

где )'ИС1 - вектор наблюдений, построенный по истинной модели (12) при полном отсутствии случайных ошибок.

Оценивание отклика проводилось по модели (1-2) при р-1:

к

У{ = е0 + еЛ + х Ру (' - ь])+ + б,- . (15)

3 21

Рассматривались модели с числом базисных функций к от 1 до 10. Для каждого значения уровня шума р проводилась серия генераций наборов данных по модели (12). Далее по моделям (15) с количествами базисных функций к = 1,2,...,10 проводилось оценивание значений отклика с использованием различных алгоритмов. По результатам оценивания определялись значения сумм квадратов (13) и (14). После окончания серии экспериментов значения сумм квадратов, соответствующих каждой комбинации «число узлов - алгоритм оценивания» усреднялись. Для удобства дальнейшего изложения введем ряд

обозначений для усредненных сумм квадратов: MS\, MS 2, MS3 - усредненные по серии наблюдений суммы квадратов (13), полученные при использовании оценок параметров (4) и алгоритмов 2, 3, 4 выбора координат узловых точек соответственно; MS А, MS5, MS6, MS 1 - усредненные по серии наблюдений суммы квадратов (13), полученные при использовании оценок параметров (5) и алгоритмов 1, 2, 3, 4 соответственно. Аналогичным образом обозначим усредненные суммы квадратов (14): MST], MST2, MST3, MST4, MST5, MST в , MST1.

В ходе вычислительных экспериментов генерация исходных данных проводилась при дисперсиях ошибок, соответствующих уровням шума р от 5 % до 25 %. Усреднение результатов проводилось по сериям из 300 наборов данных. Итоги этих экспериментов позволили сделать вывод, что при малом уровне шума (рис. 2) точность алгоритмов как с точки зрения воспроизведения исходных наблюдений, так и с точки зрения соответствия истинной модели меняется одинаковым образом. В частности, для всех алгоритмов справедливо повышение точности с увеличением количества базисных функций, при этом в среднем наилучшие результаты наблюдаются у алгоритмов, основанных на использовании обычных МНК-оценок параметров (4). Введение параметра сглаживания при малом уровне шума не влечет за собой улучшения качества модели.

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

-•-MST1 —•—MST2 —*—MST3 -MST 4 -»-MST5 -•-MST6 —1—MST7

\ \ \

\ \ \

Т-1 lV~th < 9 Ш ч к

О 2 4 6 8 10

Рис. 2. Усредненные суммы квадратов при уровне шума р = 5 %

Однако с точки зрения соответствия истинной зависимости, качество этих моделей с ростом числа базисных функций постепенно ухудшается. С этой точки зрения, при достаточном количестве базисных функций (в данном случае при к > 5), наилучшие результаты демонстрируют модели, основанные на использовании оценок (5).

Следует обратить особое внимание на модель, основанную на оценках (4) и алгоритме выбора координат узловых точек номер 2 (на графиках этой модели соответствует номер 1). Данная модель с точки зрения точности занимает промежуточное значение между моделями, использующими и не использующими сглаживание непараметрической части. С позиции точности воспроизведения исходных данных эта модель уступает моделям с номерами 2 и 3, однако в среднем оказывается точнее всех моделей, основанных на оценках (5). С позиции точности соответствия истинной зависимости в среднем эта модель хоть и уступает моделям, основанным на оценках (5), но показывает значительно лучшие результаты, чем все остальные модели, основанные на оценках (4).

—•—МЭИ -»-М5Т2 —*—М8ТЗ

-#-МБТ6 +М8Т7

.4 — -

к

Рис. 3. Усредненные суммы квадратов при уровне шума р = 20 %

Для анализа устойчивости рассматриваемых методов построения линии регрессии был проведен ряд вычислительных экспериментов с искусственными нарушениями предположений модели (12). В частности, в ходе генерации случайной ошибки в выборку включались аномальные наблюдения - «выбросы». Доля выбросов в исходных данных составляла 3 %. Роль выбросов исполняли нормально распределенные случайные величины с десятикратно

увеличенной дисперсией (1()о~ ). Координаты точек, содержащих выбросы, при каждой генерации исходных данных определялись случайно. Усредненные результаты испытаний представлены на рис. 4 и 5.

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

Рис. 4. Усредненные суммы квадратов при уровне шума р = 5 % с выбросами

Среди всех моделей, использующих сглаживание непараметрической части, наиболее устойчивой оказалась модель, в которой координаты узловых точек базисных функций определялись по второму алгоритму. Наименее устойчивыми к выбросам оказались модели, полученные при использовании алгоритмов 3 и 4. Модели, построенные с использованием алгоритма 1, показывают наилучшие результаты при увеличении числа базисных функций.

Рис. 5. Усредненные суммы квадратов при уровне шума р = 20 % с выбросами

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

ЗАКЛЮЧЕНИЕ

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

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

СПИСОК ЛИТЕРАТУРЫ

[1] Horowitz J.L. Semiparametric andNonparametric Methods in Econometrics / J.L. Horowitz. -New York: Springer, 2009. - 286 p.

[2] Ichimura H. Implementing nonparametric and semiparametric estimators Handbook ofEconometrics / H. Ichimura, P.E. Todd. - Vol. 6.-PartB.-Elsevier Science. - 2007. - P. 5369-5468.

[3] Ruppert D. Semiparametric Regression / D. Ruppert, M.P. Wand, R.J. Carroll. - New York: Cambridge university

press, 2003.-404 p.

[4] Денисов В.И. Штрафные сплайны в задаче идентификации полупараметрической регрессии / В.И. Денисов, B.C. Тимофеев, О.И. Бузмакова // Научн. вестник НГТУ. - Новосибирск: Изд-во СО РАН. - 2011. - № 4 (45). - С. 11-24.

[5] Денисов В.И. К вопросу выбора оптимальных координат узловых точек в моделях полупараметрической регрессии / В.И.Денисов, А.В. Фадцеенков // Научн. вестник НГТУ. - Новосибирск: Изд-во СО РАН. - 2012. -№4 (49).-С. 3-11.

[6] Ивахненко А.Г. Помехоустойчивость моделирования / A.T. Ивахненко, B.C. Степашко. - Киев: Наукова думка, 1985.-216 с.

[7] Кендалл М. Статистические выводы и связи / М. Кендалл, А. Стьюарт. - М.: Наука, 1973. - 466 с.

Денисов Владимир Иванович, заслуженный деятель науки РФ, доктор технических наук, профессор, академик МАН ВШ, член-корреспондент АИН РФ, профессор кафедры прикладной математики факультета прикладной математики и информатики НГТУ. Основное направление научных исследований: разработка и исследование статистических методов анализа, планирования экспериментов и прогнозирования многофакторных статистических и динамических объектов. Имеет более 250 публикаций. E-mail: videnis@nstu.ru

Тимофеев Владимир Семенович, доктор технических наук, доцент, декан факультета прикладной математики и информатики НГТУ. Основное направление научных исследований: разработка и исследование устойчивых методов и алгоритмов анализа многофакторных объектов, в том числе с использованием непараметрической статистики. Имеет более 75 публикаций, в том числе один учебник. E-mail: netsc@fpm.ami.nstu.ru

Фаддеенков Андрей Владимирович, кандидат технических наук, доцент кафедры «Теория рынка» НГТУ. Основное направление научных исследований: разработка и исследование методов и алгоритмов анализа многофакторных объектов со структурированной ошибкой. Имеет более 30 публикаций, в том числе один учебник. E-mail: fadd@fb.nstu.ru

Denisov V.I., Timofeev V.S., Faddeenkov A.V.

Investigation of algorithms of selection of knots' optimal coordinates for semiparametric models -with penalized

splines

Modifications of the regression construction methods on the basis of the semi-parametric models are suggested. New algorithms choosing of knots' optimal coordinates on the basis of criteria of the accuracy and descriptiveness of individual observations are developed. The results of comparative studies of the developed algorithms for different variants ofdata contamination taken with the computational experiments presented.

Key words: parametric and non-parametric methods, semiparametric regression, penalized splines models, basis functions, least squares method, design of experiments.