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

Научный вестник НГТУ. - 2010. - № 3(40)

УДК 62-83: 531.3

Ядерные оценки функции плотности при идентификации уравнений регрессии*

В.С. ТИМОФЕЕВ

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

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

ВВЕДЕНИЕ

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

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

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

* Получена 21 мая 2010 г.

Работа выполнена в рамках проекта № П263 ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 гг.

мального правдоподобия на основе ядерных оценок. Это позволит полностью отойти от параметрического представления распределения с целью максимально возможной «подстройки» под практически реализуемые ситуации. Первая попытка была основана на использовании характеристической функции [9]. Полученные автором результаты подтверждают хорошие адаптационные возможности алгоритма.

где X =

1. ПОСТАНОВКА ЗАДАЧИ И ОСНОВНЫЕ ПРЕДПОЛОЖЕНИЯ Рассмотрим регрессионное уравнение вида

у = х е+в,

/ 1(хи) ••• / р(х1рУ

- матрица значений регрессионных функций, имеющая пол-

(1)

/1( "^ш) / р (хщ)

ный столбцовый ранг, т. е. X) = р, е = (е^..., ер) - вектор неизвестных параметров, подлежащих оцениванию; р - количество неизвестных параметров; N - количество проведенных экспериментов; Л (х) - известные действительные функции вещественного аргумента х, Ху - заданные значения входных факторов в N наблюдениях; у = (У1,..., )Т- вектор

значений отклика; в = (в1,..., в N )т - вектор ошибок наблюдений.

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

но, что

Е (вг ) = 0, D (вг ) = а2.

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

2. МАКСИМАЛЬНОЕ ПРАВДОПОДОБИЕ И ЯДЕРНЫЕ ОЦЕНКИ ФУНКЦИИ ПЛОТНОСТИ

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

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

^жМ"?}, <2>

где X - ширина ядра (окна сглаживания), К (г) - функция ядра.

Отметим, что ядерная функция и ширина окна сглаживания должны быть выбраны таким образом, чтобы удовлетворять следующим условиям регулярности [12].

1. Четность, т. е. К (г ) = К (—г) . 1. Нт Х = 0.

N ^да

2. 0 < К (г) <да. 2. lim Ж = да.

N ^да

да

3. | К(г)dг = 1.

—да да

4. | г2К (г)dг = 1.

—да да

5. | гНК(

г) Ог <да .

Анализ литературных источников показывает [12], что существует достаточно большой выбор ядерных функций, приведем несколько наиболее популярных. 1. Квадратическая ядерная функция

К (г ) =— (1 — г2 ), У| г| < 1 и К (г) = 0 иначе.

2. Ядро Епанечникова

К(г) = 3(1 — г2) , У|г| < 1 и К(г) = 0 иначе.

3. Ядро Гаусса

1 ( —г 2 ^

К(г) = "/1ГехР —Г . V2* ^ 2 )

4. Треугольное ядро

К (г ) = (1 — |г|) , У| г| < 1 и К (г) = 0 иначе.

5. Прямоугольное ядро

К(г) = -1, У|г < 1 и К (г ) = 0 иначе.

6. Квадратичная ядерная функция

3 3

К(г) = —=--= г2, У|г\ <75 и К(г) = 0 иначе.

4 7 4л/5 20л/5 11

7. Ядро, предложенное С.А. Айвазяном в [1],

К ( г )= 1

1 + г2 '

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

—оо

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

Далее будем считать, что выбор X обеспечивает унимодальность распределения. Кроме того, в силу предположения о независимости случайных ошибок уравнения (1) значения остатков еI = уI — х10 (х^ - / -я строка матрицы X из (1)) также будут статистически независимыми случайными величинами с плотностью распределения у (и, 0). Тогда для оценивания

параметров уравнения (1) можно воспользоваться методом максимального правдоподобия [5, 6, 8]. Учитывая, что остатки наблюдаемы, т. е. их значения определяются на основе имеющихся исходных данных, можно записать логарифмическую функцию правдоподобия

/ ( К I N I I

I(еь...,ем,0) = 1п е,0) =11п(у(е,0)) . (3)

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

3. АЛГОРИТМ ОЦЕНИВАНИЯ ПАРАМЕТРОВ

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

Шаг 1. Определение начального приближения для вектора неизвестных параметров уравнения (1), в качестве которого автором была использована оценка метода наименьших квадратов [5], что позволило сократить число итераций и время вычислений.

Шаг 2. Вычисление остатков регрессионного уравнения.

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

—1

Х = N 5 вычисление значений искомой функции плотности у(и1) по соотношению (2) в точках, соответствующих вычисленным значениям остатков е^ .

Шаг 4. Вычисление значения логарифмической функции правдоподобия (3).

Шаг 5. Поиск очередного приближения оценок вектора неизвестных параметров 0

0к+1 = а^ тах I(е1,е2,...,eN,0к).

Шаг 6. Если

0 к+1 — к

< в , то завершение процесса, в противном случае к := к +1 и

переход на шаг 2 (в - заданная погрешность вычисления).

4. РЕЗУЛЬТАТЫ ВЫЧИСЛИТЕЛЬНЫХ ЭКСПЕРИМЕНТОВ

Для исследования разработанного алгоритма оценивания вектора неизвестных параметров 0 уравнения (1) проводились многочисленные вычислительные эксперименты. Приведем лишь некоторые из полученных результатов. В качестве истинной зависимости рассмотрим уравнение регрессии

у = 0О +01х + 02 х2 + в, (4)

где количество регрессоров р = 3, значения входных факторов х ^ выбирались из отрезка [—1,1], истинные значения неизвестных параметров 90 = 50, 01 = 25, 02 = 10. Случайные ошибки в 7 моделировались независимыми и одинаково распределенными с функцией распределения вида

F(х) = (1 — X)F1 (х, т1, с 1 ) + А^2 (х, т2, с 2), (5)

где F7 (х, т7, с 7) - функция нормального распределения с математическим ожиданием, равным т7, и дисперсией с2 , 7 = 1,2, X е[0,1] - параметр смеси. Во всех проведенных вычислительных экспериментах т^ = т2 = 0 (если по тексту не оговорено противное).

Представление (5) позволяет моделировать ошибку с различной степенью отклонения от нормального распределения, в том числе появление отдельных, довольно грубых засоряющих наблюдений - «выбросов». Параметр X определяет соответствующие доли наблюдений с 22

дисперсиями С1 и С2 в выборке. При X = 0 и X = 1 ошибка будет иметь нормальное распределение. В проведенных вычислительных экспериментах полагалось, что с2 > С. Однако при

22

моделировании задавались не сами значения дисперсий С1 и С2 , а им соответствующие значения уровня шума. Уровень шума введен в [7] и определяется как отношение «шум»/«сигнал» в %:

р = -100%, с

1 "

где с - дисперсия ошибки; с =-^ (у° — у ) - интенсивность сигнала (не зашумлен-

П — 17=1

ных измерений у0 ).

В качестве показателей точности оценивания параметров использовались Ll нормы отклонений оценок неизвестных параметров от истинных значений

£ =

0 ист — 0

0

ист

= 0ист — 0

Прежде всего проиллюстрируем качество восстановления функции плотности при разном величине ширины окна X. На рис. 1 представлен результат восстановления функции плотности с помощью соотношения (2) для нормально распределенной ошибки при N = 276 и X = 0.3,1.0,1.8.

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

Далее рассмотрим результаты оценивания параметров тестового регрессионного уравнения (4). Для различных комбинаций X и р проводилось по 600 вычислительных экспериментов. Каждый такой эксперимент заключался в моделировании выборки исходных данных в соответствии с моделью (4) с последующим оцениванием параметров этой модели предлагаемым в этой статье алгоритмом, алгоритмом, основанным на восстановлении функции плотности через характеристическую функцию [9], а также методом наименьших квадратов [5] и

и

знаковым методом [2, 4]. В качестве итоговых показателей точности оценивания использовались усредненные по 600 экспериментам значения показателей ^ и 52 .

а б

в

Рис. 1. Функции плотности остатков: при а - Х = 0.3 б - Х = 1.0 в - Х = 1.8

В таблице представлены результаты, полученные для нормально распределенной ошибки наблюдения (X = 0) при разном объеме выборки. Дисперсия ошибки соответствовала уровню шума 5 %.

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

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

Метод оценивания Объем выборки, N

100 200 500 100 200 500

51 52

МНК 0,0187 0,021 0,013 0,404 0,452 0,293

Ядерная 0,025 0,017 0,010 0,538 0,368 0,228

Характеристическая функция 0,025 0,018 0,011 0,527 0,434 0,297

Знаковый метод 0,019 0,013 0,008 0,411 0,301 0,187

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

Далее рассмотрим изменение точности оценивания в зависимости от дисперсии (уровня шума) нормально распределенной ошибки (X = 0). Для этого будем последовательно изменять

2

уровень шума, соответствующий дисперсии с , от 2.5 до 40 % с шагом 2.5. Результаты проведенных экспериментов приведены на рис. 2.

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

Рис. 2. Точность оценивания в зависимости от уровня шума ( N = 100 )

Также проведено исследование точности оценивания вектора неизвестных параметров 9 при разной степени отклонения распределения случайной ошибки от нормального распределения. Для этого изменению подвергался параметр смеси X . При малых значениях X в выборке будет появляться небольшое число выбросов, а при значениях X, близких к 0.5, можно говорить об изменении формы распределения. Было зафиксировано Р1 = 5%, Р2 = 50%, а доля

выбросов X изменялась от 0 до 0.5 с шагом 0.02. Результаты оценивания представлены на рис. 3, где показано изменение показателя для предложенного алгоритма, алгоритма, основанного на восстановлении функции плотности через характеристическую функцию, знакового метода и МНК. Объем выборки 200.

Рис. 3. Точность оценивания в зависимости от X ( N = 200 )

Из рис. 3 видно, что метод наименьших квадратов заметно уступает в точности оценивания всем остальным методам, поэтому наиболее интересно провести анализ точности оценивания оставшихся методов. Разработанный здесь алгоритм при всех рассмотренных уровнях засорения уверенно превосходит по точности алгоритм, основанный на восстановлении функции плотности через характеристическую функцию. Кроме того, на средних уровнях засорения предложенный алгоритм показывает примерно одинаковую со знаковым методом точность оценивания. Как известно, знаковый метод относится к робастным методам и имеет хорошую устойчивость к наличию в выборке определенной доли грубых ошибок, что показано в [2, 4]. При этом на малой и очень большой степенях засорения знаковый метод несколько точнее. Последнее можно объяснить потерей свойства унимодальности восстановленной плотности, ведь именно малые и большие степени засорения наиболее подвержены к появлению «дополнительных» максимумов у функции плотности.

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

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

также был проведен ряд вычислительных экспериментов, в которых зафиксировано Р1 = 5 %, Р2 = 50 %, доля выбросов X изменялась от 0 до 0.5 с шагом 0.02. Объем выборки 100, 200, 500. Полученные результаты представлены на рис. 4, где показано изменение отношения показателя для знакового метода к значению этого же показателя, но рассчитанного для пред-

Рис. 4. Сравнение точности оценивания при разных объемах выборки

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

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

ЗАКЛЮЧЕНИЕ

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

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

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

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

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

[1] Айвазян С.А. Прикладная статистика. Исследование зависимостей: справ. изд. - М.: Финансы и статистика,

1985.

[2] Болдин М.В., Симонова Г.И., Тюрин Ю.Н. Знаковый статистический анализ линейных моделей. - М.: Наука Физматлит, 1997.

[3] Денисов В.И., Тимофеев В.С. Оценивание параметров регрессионных зависимостей с использованием аппроксимации Грамма-Шарлье // Автометрия. - 2008.- Т. 44, № 6. - С. 3-12.

[4] Денисов В.И., Тимофеев В.С. Знаковый метод: преимущества, проблемы, алгоритмы // Науч. вестн. НГТУ. - 2001. - №1(10). - С. 21-35.

[5] Дрейпер Н., Смит Н. Прикладной регрессионный анализ. - М.: Статистика, 1973.

[6] Закс Ш. Теория статистических выводов. - М.: Мир, 1975.

[7] Ивахненко А.Г., Степашко В.С. Помехоустойчивость моделирования. - Киев: Наук. думка, 1985.

[8] Кендалл М., Стьарт А. Теория распределений. - М.: Наука, 1966.

[9] Тимофеев В.С. Оценивание параметров регрессионных зависимостей на основе характеристической функции // Науч. вестн. НГТУ. - 2010. - № 2(39). - С. 43-52.

[10] Тимофеев В.С. Оценивание параметров регрессионных зависимостей с использованием кривых Пирсона. Ч. 1 // Науч. вестн. НГТУ. - 2009. - № 4(37). - С. 57-66.

[11] Химмельблау Д. Прикладное нелинейное программирование. - М.: 1975.

[12] Pagan A., Ullah A. Nonparametric econometrics. - N. Y., 1999.

Тимофеев Владимир Семенович, кандидат технических наук, доцент кафедры теории рынка Новосибирского государственного технического университета. Основное направление научных исследований - разработка и исследование устойчивых методов и алгоритмов анализа многофакторных объектов, в том числе с использованием непараметрической статистики. Имеет более 50 публикаций, в том числе один учебник.

E-mail: [email protected] Тел. 8(383) 346-31-72

V.S. Timofeev

The kernel estimation of density function in the regression identification problem

The problem of the parameters estimation for regression models is investigated. The probability density function in maximum likelihood method reconstructed by use a kernel functions is proposed. The advantages of the proposed approach are notified. The results of the computing experiments are discussed. The precision of the method is compare with other algorithms.

Key words: regression equation, parameters estimation, kernel function, maximum likelihood method, density function.