Алгоритм LOWESS

Материал из MachineLearning.

(Различия между версиями)
Перейти к: навигация, поиск
(Литература)
(Введение)
 
(181 промежуточная версия не показана)
Строка 1: Строка 1:
-
{{stop|
+
== Введение ==
-
'''Статья плохо доработана.'''<br/>
+
-
Имеются указания по её улучшению:<br/>
+
-
[[Изображение:Edit1.jpg|||100px]][[Изображение:Edit2.jpg|||100px]]
+
-
}}
+
-
'''Алгоритм LOWESS''' ''(locally weighted scatter plot smoothing)'' - локально взвешенное сглаживание.
+
-
=== Постановка задачи ===
+
[[Изображение:Loess_smooth.jpg|thumb|right|300px|Рис. 1. Пример применения lowess-сглаживания]]
-
:Решается задача восстановления регрессии. Задано пространство объектов <tex>X</tex> и множество возможных
+
-
ответов <tex>Y=R</tex>. Существует неизвестная целевая зависимость <tex> y^*: X \rightarrow Y</tex>,
+
-
значения которой известны только на объектах обучающей выборки <tex> X^m={(x_i, y_i)}_{i=1}^m</tex>.
+
-
Требуется построить алгоритм <tex>a: X \rightarrow Y </tex>, аппроксимирующий целевую зависимость <tex>y^*</tex>.
+
-
=== Непараметрическая регрессия ===
+
: Данная методика была предложена Кливлендом(Cleveland) в 1979 году для моделирования и сглаживания двумерных данных <tex> X^m={(x_i, y_i)}_{i=1}^m</tex>. Эта техника предоставляет общий и гибкий подход для приближения двумерных данных.
-
:Непараметрическое восстановление регрессии основано на идее, что значение <tex>a(x)</tex> вычисляется
+
-
для каждого объекта <tex>x</tex> по нескольким ближайшим к нему объектам обучающей выборки.
+
-
В формуле Надарая–Ватсона для учета близости объектов <tex>x_i</tex> обучающей выборки к объекту <tex>a(x)</tex>
+
: '''Локально-линейная модель''' loess(lowess) может быть записана в виде:
-
предлагалось использовать невозрастающую, гладкую, ограниченную функцию
+
:: <tex> y_t=\alpha_t+\beta_t x_t + \varepsilon_t.</tex>
-
<tex>K: [0,\infty) \rightarrow [0,\infty)</tex>, называемую ядром:
+
-
::<tex>w_i(x) = K\left( \frac{\rho(x, x_i)}{h}\right)</tex>
+
: Эта модель может быть '''расширена''' на случай '''локально-квадратичной зависимости''' и на модель с '''бо‘льшим числом независимых переменных'''.
-
Параметр <tex>h</tex> называется ''шириной ядра'' или ''шириной окна сглаживания''. Чем меньше <tex>h</tex>,
+
: Параметры <tex>\alpha_t</tex> и <tex>\beta_t</tex> локально линейной модели оцениваются с помощью локально взвешенной регрессии, которая присваивает объекту тем больший вес, чем более близок он к объекту <tex>t</tex>.
-
тем быстрее будут убывать веса <tex>w_i(x)</tex> по мере удаления <tex>x_i</tex> от <tex>x</tex>.
+
-
В общем случае <tex>h</tex> зависит от объекта <tex>x</tex>, т.е. <tex>h=h(x)</tex>. Тогда веса вычисляются по формуле
+
-
<tex>\textstyle w_i(x) = K\left( \frac{\rho(x, x_i)}{h(x)}\right)</tex>
+
-
====Оптимизация ширины окна====
+
:Степень сглаживания определяется '''параметром сглаживания''' <tex>f</tex>, который выбирает пользователь.
-
Чтобы оценить при данном <tex>h</tex> и <tex>K</tex> точность локальной аппроксимации в точке <tex>x_i</tex>,
+
-
саму эту точку необходимо исключить из обучающей выборки. Если этого не делать, минимум ошибки будет
+
-
достигаться при <tex>h\rightarrow 0</tex>. Такой способ оценивания оптимальной ширины окна называется скользящим контролем
+
-
''с исключением объектов по одному'' (leave-one-out, LOO):
+
-
::<tex>LOO(h,X^m) = \sum_{i=1}^m{\left(a_h(x_i;X^m\setminus\{x_i\}) - y_i \right)^2} \rightarrow min\limits_h</tex>
+
: Параметр <tex>f</tex> '''указывает, какая доля (fraction) данных''' используется в процедуре. Если <tex>f = 0.5</tex>, то только половина данных используется для оценки и влияет на результат, и тогда мы получим умеренное сглаживание. С другой стороны, если <tex>f = 0.8</tex>, то используются восемьдесят процентов данных, и сглаживание намного сильнее. Во всех случаях веса данных тем больше, чем они ближе к объекту <tex>t</tex>.
 +
: Процедура оценки использует '''не''' метод наименьших квадратов, а более устойчивый ( робастный ) метод, который принимает меры против выбросов.
 +
:График приближенных значений
 +
:: <tex> \hat{y_t}=\hat{\alpha_t}+\hat{\beta_t}x_t </tex>
 +
:от <tex>x_t</tex> '''полезен для принятия решения о характере связи между <tex>y_t</tex> и <tex>x_t</tex>'''. Для проверки качества приближения полученного с помощью процедуры устойчивого loess полезно посмотреть на график остатков обычной регрессии, то есть в осях '''(i)''' остатки от числа наблюдения '''(ii)''' остатки от прибли‘женных значений, '''(iii)''' остатки от значений независимой переменной. Как показал Кливленд, может быть предпочтительно использовать график в осях модули остатков от полученных приближенных значений вместо графика '''(ii)''' для устойчивого loess сглаживания, чтобы проверить наличие тренда или других систематических особенностей.
-
== Проблема выбросов ==
+
:Когда <tex>m > 100</tex> вычисления могут быть слишком долгими, в этом случае можно сократить количество вычислений, оценивая <tex>\hat{\alpha_t}</tex> и <tex>\hat{\beta_t}</tex> только в точках отстоящих друг от друга как минимум на <tex>\delta</tex> единиц, где параметр <tex>\delta</tex> может задаваться либо приниматься по умолчанию. Рекомендуемые значения
-
:Оценка Надарайя–Ватсона <tex>\textstyle a_h(x,X^m) = \frac{\sum_{i=1}^m{y_iw_i}}{\sum_{i=1}^m{w_i}}</tex>
+
:: <tex>\delta=0, </tex> если <tex>m <= 100 </tex>
-
'''крайне чувствительна к большим одиночным выбросам'''. На практике легко идентифицируются только грубые ошибки,
+
:: <tex>\delta=0.03*IQR,</tex> если <tex>m > 100</tex>, где <tex>IQR</tex> — ['''межквартильный размах'''](Interquartile range).
-
возникающие, например, в результате сбоя оборудования или невнимательности персонала при подготовке данных.
+
:С такими параметрами вычисления будут выполнены для примерно 100 точек.
-
В общем случае можно лишь утверждать, что чем больше величина ошибки
+
-
::<tex>\varepsilon_i = \left | a_h \left (x_i;X^m\setminus\{x_i\} \right) -y_i \right |</tex>
+
=== Примеры ===
 +
[[Изображение:Loess_f_s.jpg|thumb|Рис. 2. Задание параметра сглаживания <tex>f</tex>|300px]]
 +
:На '''Рис. 2'''. Приведена иллюстрация уровня сглаживания в зависимости от значения параметра <tex>f</tex>
-
тем в большей степени прецедент <tex>(x_i,y_i)</tex> является выбросом , и тем меньше должен быть его вес.
+
:Сглаживание также может быть локально квадратичным, в этом случае модель для <tex>y_t</tex> имеет вид
-
Эти соображения приводят к идее домножить веса <tex>w_i(x)</tex> на коэффициенты
+
:: <tex> y_t=\alpha_t+\beta_t x_t +\gamma_t x_t^2+ \varepsilon_t.</tex>
-
<tex>\gamma_i = \bar{K}(\varepsilon_i)</tex>, где <tex>\bar{K}</tex> — ещё одно ядро, вообще говоря,
+
Примеры сглаживания с '''квадратичным локальным приближением''' показаны на '''Рис. 3'''.
-
отличное от <tex>K</tex>.
+
 
 +
[[Изображение:Loess_quadr.jpg|thumb|Рис. 3. Локально квадратичное сглаживание|300px]]
 +
 
 +
== Технические детали алгоритма ==
 +
Базовое предположение состоит в следующем
 +
::<tex>y_t=g(x_t)+\varepsilon_t , t=1,\ldots,m</tex>
 +
где <tex>g(x)</tex> - функция сглаживания, остатки <tex>\varepsilon_t</tex> имеют нулевое математическое ожидание и фиксированную дисперсию. Затем сглаживание <tex>g</tex> мы приближаем локально-линейной (локально квадратичной, в случае нелинейной модели) функцией, чтобы получить
 +
:: <tex> y_t=\alpha_t + \beta_t x_t + \varepsilon_t</tex>.
 +
Для четкого определения алгоритма поясним концепцию '''локальных весов''' <tex>w(x_t)</tex> и '''робастных весов''' <tex>\delta(x_t)</tex>.
 +
=== Локальные веса ===
 +
:Рассмотрим один из широко распространенных примеров – функцию
 +
::<tex>W(z)=(1-|z|^3)^3, \, \, |z|<=1 \\ W(z)=0, \,\, |z|>1 </tex>
 +
:Для заданного параметра <tex>0 < f < 1</tex> пусть <tex>r</tex> - ближайшее целое число к произведению <tex> f*m</tex>. Пусть <tex> h_t</tex> расстояние до <tex>r</tex>-того ближайшего соседа объекта <tex>x_t</tex>. Тогда локальный вес для любого объекта <tex>x</tex> в окрестности <tex>x_t</tex> есть
 +
::<tex>w(x)=W\left(\frac{x-x_t}{h_t}\right)</tex>.
 +
==== Замечание ====
 +
:Более общий подход к определению локальных весов состоит в выборе ширины окна <tex>h</tex>, в общем случае <tex>h=h(x_t)</tex>, то есть зависящей от объекта <tex>x_t</tex>, и ядровой функции <tex>K(x)=K\left(\frac{\rho(x,x_t)}{h(x_t)} \right )</tex>. Тогда локальные веса вычисляются по формулам
 +
::<tex>w(x)=K \left( \frac{\rho(x,x_t)}{h(x_t)} \right ).</tex>
 +
:В этом случае отпадает необходимость задания параметра сглаживания <tex>f</tex> и его смысл эквивалентен выбору ширины окна <tex>h=h(x_t)</tex>.
 +
 
 +
=== Робастные веса ===
 +
Пусть
 +
::<tex>X^m\setminus\{x_t\}</tex> – обучающая выборка за исключением элемента <tex>x_i</tex>,
 +
::<tex>\hat{y_t}:=a \( x_t;X^m\setminus\{x_t\} \)</tex> – ответ алгоритма <tex>a</tex>, обученного на выборке <tex>X^m\setminus\{x_t\}</tex> при работе на объекте <tex>x_t</tex>.
 +
::<tex>\hat{\varepsilon_t}= \|\hat{y_t}-y_t \|</tex> – ошибка алгоритма на объекте <tex>x_i</tex>('''ошибка скользящего контроля''').
 +
Пусть <tex>s </tex> - есть медиана величин<tex>\hat{\varepsilon_1},\ldots,\hat{\varepsilon_m}</tex>.
 +
тогда <tex>\delta_t = \bar{K}(\frac{\varepsilon_t}{6s})</tex>, где
 +
::<tex>\bar{K}(z)=\left{ (1-|z|^2)^2, \, \, \, |z|\le 1 \\ 0, \, \, \, |z|>1 \right. </tex>
 +
 
 +
==== Замечание ====
 +
:Возможны и другие варианты выбора весов <tex>\delta_t</tex>, например, занулить <tex>p</tex> штук, соответствующих наибольшим <tex>\hat{\varepsilon_t}</tex>. Это соответствует ядру
 +
::<tex>\bar{K}(z)=[z<\hat{\varepsilon}^{(m-p)}],</tex>
 +
где <tex>\hat{\varepsilon}^{(m-p)}</tex> –- <tex>(m-p)</tex> - тый член вариационного ряда <tex>\hat{\varepsilon}^{(1)}<=\ldots<=\hat{\varepsilon}^{(m)}</tex>
 +
:В качестве весовой ядерной функции можно взять функцию Хубера (Huber, 1964) на которой основаны *[M-оценки]
 +
<tex>K(z)= \left{ z^2, \, \, \, |z| \le c\\ 2c|z|-c_2, \, \, \, |z|>c \right. </tex>
 +
Чтобы вычислить <tex>K(z)</tex> необходимо выбрать параметр устойчивости <tex>c</tex>. Одно популярное прикладное правило – <tex>c = 1,345 * s</tex> , где <tex>s</tex> – '''робастная мера масштаба''', такая как медианное абсолютное отклонение от медианы (MAD). Это популярное правило обеспечивает 95%-ую эффективность относительно гомоскедастичной нормальной модели в проблеме местоположения.
=== Алгоритм LOWESS ===
=== Алгоритм LOWESS ===
==== Вход ====
==== Вход ====
-
<tex>X^m</tex> - обучающая выборка;
+
:<tex>X^m</tex> - обучающая выборка;
-
<tex>w_i \, i=1,\ldots,m</tex> весовые функции;
+
:<tex>w_t, \,\,\, t=1,\ldots,m</tex> весовые функции;
==== Выход ====
==== Выход ====
-
Коэффициенты <tex>\gamma_i, i=1,\ldots,m</tex>
+
Коэффициенты <tex>\delta_t, \,\,\, t=1,\ldots,m</tex>
-
==== Алгоритм ====
+
==== Алгоритм 1.1====
-
:1: инициализация
+
:'''1.''' Построить линейную регрессию во всех <tex>t=1,\ldots,t=m</tex> точках, используя весовые функции <tex>w_t</tex>, тем самым получим оценки для параметров модели <tex>\hat{\alpha_t}, \hat{\beta_t}</tex>.
-
::<tex>\gamma_i:=1, i=1,\ldots,m</tex>
+
:А также приближения <tex>\hat{y_t}=\hat{\alpha_t}+\hat{\beta_t}x_t</tex>.
-
:2: '''повторять'''
+
:'''2.''' Инициализируем остатки <tex>\hat{\varepsilon_t}= \| \hat{y_t} - y_t \|</tex>. Вычислим робастные веса <tex>\delta_t =\bar{K}(\hat{\varepsilon_t})</tex>
-
:3: вычислить оценки скользящего контроля на каждом объекте:
+
:'''3.''' '''повторять'''
-
::<tex>a_i:=a_h\( x_i;X^m\setminus\{x_i\} \)=\frac{\sum_{j=1,j\ne i}^m y_j\gamma_j w_j}{\sum_{j=1,j\ne i}^m \gamma_j w_j },\;i=1,\ldots,m</tex>
+
::'''4.''' Построить линейную регрессию во всех <tex>t=1,\ldots,t=m</tex> точках, используя весовые функции <tex>\delta_t w_t</tex>, тем самым получим оценки для параметров модели <tex>\hat{\alpha_t}, \hat{\beta_t}</tex>. А также приближения <tex>\hat{y_t}=\hat{\alpha_t}+\hat{\beta_t}x_t</tex>.
-
:4: вычислить новые значения коэффициентов <tex>\gamma_i</tex>:
+
::'''5.''' По новому набору значений <tex>\hat{\varepsilon_t}= \| \hat{y_t} - y_t \|</tex> вычислить новые значения коэффициентов <tex>\delta_t =\bar{K}(\hat{\varepsilon_t})</tex>.
-
::<tex>\varepsilon_i = \left | a_i -y_i \right |</tex>
+
:'''6.''' '''пока''' веса <tex>\delta_t</tex> не стабилизируются
-
::<tex>\gamma_i:=\bar{K}( \varepsilon_i ) ,\;i=1,\ldots,m</tex>;
+
-
:5: '''пока''' коэффициенты <tex>\gamma_i</tex> не стабилизируются
+
-
Коэффициенты <tex>\gamma_i</tex>, как и ошибки <tex>\varepsilon_i</tex>, зависят от функции <tex>a_h</tex>, которая,
+
При использовании ядровых функций для оценки локальных весов объектов и робастных весов алгоритм модифицируется следующим образом:
-
в свою очередь, зависит от <tex>\gamma_i</tex>. На каждой итерации строится функция <tex>a_h</tex>,
+
-
затем уточняются весовые множители <tex>\gamma_i</tex>. Как правило, этот процесс сходится довольно быстро.
+
-
Он называется локально взвешенным сглаживанием (locally weighted scatter plot smoothing, LOWESS).
+
-
=== Выбор ядра <tex>\bar{K} </tex>===
+
==== Алгоритм 1.2====
 +
:'''1.''' Инициализировать <tex>\delta_1:=\ldots=\delta_m:=1</tex>
 +
:'''2.''' '''повторять'''
 +
::'''3.''' Вычислить оценки скользящего контроля на каждом объекте
 +
:: <tex> \hat{y_t}:=a(x_t; X\setminus\{ x_t\}) = \frac{ \sum_{i=1, i\neq t }^{m} {y_i \delta_i K\left( \frac{\rho(x_i,x_t)}{h(x_t)}\right)} } {\sum_{i=1, i\neq t }^{m} {y_i K\left( \frac{\rho(x_i,x_t)}{h(x_t)}\right)} }</tex>
 +
::'''4.''' По набору значений <tex>\hat{\varepsilon_t}= \| \hat{y_t} - y_t \|</tex> вычислить новые значения коэффициентов <tex>\delta_t</tex>.
 +
:'''5.''' '''пока''' веса <tex>\delta_t</tex> не стабилизируются
-
:В качестве ядра <tex>\bar{K} </tex> большинство практических источников рекомендуют использовать следующее:
 
-
Пусть <tex>s </tex> - есть медиана коэффициентов <tex>\gamma_1,\ldots,\gamma_m</tex>,
+
Коэффициенты <tex>\delta_t</tex>, как и ошибки <tex>\varepsilon_t</tex>, зависят от функции <tex>a</tex>, которая,
-
тогда <tex>\gamma_i = K(\frac{\varepsilon_i}{6s})</tex>, где
+
в свою очередь, зависит от <tex>\delta_t</tex>. На каждой итерации строится функция <tex>a</tex>,
-
::<tex>K(z)=(1-|z|^2)^2, \, \, |z|<=1 \\ K(z)=0, \,\, |z|>1 </tex>
+
затем уточняются весовые множители <tex>\delta_t</tex>. Как правило, этот процесс сходится довольно быстро.
 +
'''Однако в практических реализациях имеет смысл вводить ограничение на количество итераций, как правило, это 2-3 итерации'''.
-
Более простой вариант, состоит в отбросе <tex>t</tex> коэффициентов, соответствующих объектам с максимальными <tex>\varepsilon_i</tex>. Это соотвествует ядру
+
== Примеры==
-
::<tex>K(z)=[z<\varepsilon^{(m-t)}], \, \, \\ K(z)=0, \,\, |z|>=\varepsilon^{(m-t)}, </tex>
+
[[Изображение:Loess_sample.jpg|thumb|Рис. 4. Пример для модельных данных|300px]]
-
где <tex>\varepsilon^{(m-t)}</tex> –- <tex>(m-t)</tex> - тый член вариационного ряда <tex>\varepsilon^{(1)}<=,\ldots,<=\varepsilon^{(m)}</tex>
+
:На '''рисунке 4''' представлен пример '''робастного локально-линейного сглаживания''' с помощь алгоритма LOWESS. С числом итераций цикла равным 2 и параметром сглаживания <tex> f=0.5</tex>, то есть для приближения используется <tex> k=\left\lfloor n/2 \right\rfloor </tex> ближайших точек выборки.
-
== Примеры применения==
+
'''Конечно, подход вычислительно достаточно требовательный, однако этот метод заслуживает внимания тех исследователей, которые обеспокоены наличием выбросов в данных. В частности он активно применяется в биологии в области генетических исследований.'''
-
==Литература==
+
== Литература ==
-
# {{книга
+
-
|автор = Воронцов К.В.
+
-
|заглавие = Лекции по алгоритмам восстановления регрессии
+
-
|год = 2007
+
-
|ссылка = http://www.ccas.ru/voron/download/Regression.pdf
+
-
}}
+
# {{книга
# {{книга
Строка 101: Строка 116:
|ссылка = http://www.stats.uwo.ca/faculty/aim/2004/04-259/notes/RL.pdf
|ссылка = http://www.stats.uwo.ca/faculty/aim/2004/04-259/notes/RL.pdf
}}
}}
-
 
+
# {{книга
 +
|автор = Хардле В.
 +
|заглавие = Прикладная непараметрическая регрессия.
 +
|издательство = Мир
 +
|год = 1993
 +
}}
 +
# {{книга
 +
|автор = Воронцов К.В.
 +
|заглавие = Лекции по алгоритмам восстановления регрессии
 +
|год = 2007
 +
|ссылка = http://www.ccas.ru/voron/download/Regression.pdf
 +
}}
# {{книга
# {{книга
|автор = John A Berger, Sampsa Hautaniemi, Anna-Kaarina Järvinen, Henrik Edgren, Sanjit K Mitra and Jaakko Astola
|автор = John A Berger, Sampsa Hautaniemi, Anna-Kaarina Järvinen, Henrik Edgren, Sanjit K Mitra and Jaakko Astola
Строка 108: Строка 134:
|год = 2004
|год = 2004
|ссылка = http://www.biomedcentral.com/1471-2105/5/194
|ссылка = http://www.biomedcentral.com/1471-2105/5/194
 +
}}
 +
# {{книга
 +
|автор = Maronna, A., R. Martin, V. Yohai
 +
|заглавие = Robust Statistics: Theory and Methods.
 +
|издательство = Wiley
 +
|год = 2006
 +
}}
 +
# {{книга
 +
|автор = Расин, Джеффри
 +
|заглавие = «Непараметрическая эконометрика: вводный курс»
 +
|издательство = Квантиль, №4, стр. 7–56.
 +
|год = 2008
 +
}}
 +
# {{книга
 +
|автор = Huber, P.J.
 +
|заглавие = Robust estimation of a location parameter. Annals of Statistics 35
 +
|страницы = 73–101
 +
|год = 1964
}}
}}
Строка 114: Строка 158:
* [[Регрессионный анализ]]
* [[Регрессионный анализ]]
* [http://en.wikipedia.org/wiki/Local_regression Local regression]
* [http://en.wikipedia.org/wiki/Local_regression Local regression]
-
* Расин, Джеффри (2008) «Непараметрическая эконометрика: вводный курс», Квантиль, №4, стр. 7–56.
 
-
 
-
[[Категория:Регрессионный анализ]]
 
-
{{Задание|Валентин Голодов|Vokov|31 декабря 2009}}
+
[[Категория:Непараметрическая регрессия]]
-
+
[[Категория:Робастная регрессия]]

Текущая версия

Содержание

Введение

Рис. 1. Пример применения lowess-сглаживания
Рис. 1. Пример применения lowess-сглаживания
Данная методика была предложена Кливлендом(Cleveland) в 1979 году для моделирования и сглаживания двумерных данных  X^m={(x_i, y_i)}_{i=1}^m. Эта техника предоставляет общий и гибкий подход для приближения двумерных данных.
Локально-линейная модель loess(lowess) может быть записана в виде:
 y_t=\alpha_t+\beta_t x_t + \varepsilon_t.
Эта модель может быть расширена на случай локально-квадратичной зависимости и на модель с бо‘льшим числом независимых переменных.
Параметры \alpha_t и \beta_t локально линейной модели оцениваются с помощью локально взвешенной регрессии, которая присваивает объекту тем больший вес, чем более близок он к объекту t.
Степень сглаживания определяется параметром сглаживания f, который выбирает пользователь.
Параметр f указывает, какая доля (fraction) данных используется в процедуре. Если f = 0.5, то только половина данных используется для оценки и влияет на результат, и тогда мы получим умеренное сглаживание. С другой стороны, если f = 0.8, то используются восемьдесят процентов данных, и сглаживание намного сильнее. Во всех случаях веса данных тем больше, чем они ближе к объекту t.
Процедура оценки использует не метод наименьших квадратов, а более устойчивый ( робастный ) метод, который принимает меры против выбросов.
График приближенных значений
 \hat{y_t}=\hat{\alpha_t}+\hat{\beta_t}x_t
от x_t полезен для принятия решения о характере связи между y_t и x_t. Для проверки качества приближения полученного с помощью процедуры устойчивого loess полезно посмотреть на график остатков обычной регрессии, то есть в осях (i) остатки от числа наблюдения (ii) остатки от прибли‘женных значений, (iii) остатки от значений независимой переменной. Как показал Кливленд, может быть предпочтительно использовать график в осях модули остатков от полученных приближенных значений вместо графика (ii) для устойчивого loess сглаживания, чтобы проверить наличие тренда или других систематических особенностей.
Когда m > 100 вычисления могут быть слишком долгими, в этом случае можно сократить количество вычислений, оценивая \hat{\alpha_t} и \hat{\beta_t} только в точках отстоящих друг от друга как минимум на \delta единиц, где параметр \delta может задаваться либо приниматься по умолчанию. Рекомендуемые значения
\delta=0, если m <= 100
\delta=0.03*IQR, если m > 100, где IQR — [межквартильный размах](Interquartile range).
С такими параметрами вычисления будут выполнены для примерно 100 точек.

Примеры

Рис. 2. Задание параметра сглаживания
Рис. 2. Задание параметра сглаживания f
На Рис. 2. Приведена иллюстрация уровня сглаживания в зависимости от значения параметра f
Сглаживание также может быть локально квадратичным, в этом случае модель для y_t имеет вид
 y_t=\alpha_t+\beta_t x_t +\gamma_t x_t^2+ \varepsilon_t.

Примеры сглаживания с квадратичным локальным приближением показаны на Рис. 3.

Рис. 3. Локально квадратичное сглаживание
Рис. 3. Локально квадратичное сглаживание

Технические детали алгоритма

Базовое предположение состоит в следующем

y_t=g(x_t)+\varepsilon_t , t=1,\ldots,m

где g(x) - функция сглаживания, остатки \varepsilon_t имеют нулевое математическое ожидание и фиксированную дисперсию. Затем сглаживание g мы приближаем локально-линейной (локально квадратичной, в случае нелинейной модели) функцией, чтобы получить

 y_t=\alpha_t + \beta_t x_t + \varepsilon_t.

Для четкого определения алгоритма поясним концепцию локальных весов w(x_t) и робастных весов \delta(x_t).

Локальные веса

Рассмотрим один из широко распространенных примеров – функцию
W(z)=(1-|z|^3)^3, \, \, |z|<=1 \\ W(z)=0, \,\, |z|>1
Для заданного параметра 0 < f < 1 пусть r - ближайшее целое число к произведению  f*m. Пусть  h_t расстояние до r-того ближайшего соседа объекта x_t. Тогда локальный вес для любого объекта x в окрестности x_t есть
w(x)=W\left(\frac{x-x_t}{h_t}\right).

Замечание

Более общий подход к определению локальных весов состоит в выборе ширины окна h, в общем случае h=h(x_t), то есть зависящей от объекта x_t, и ядровой функции K(x)=K\left(\frac{\rho(x,x_t)}{h(x_t)} \right ). Тогда локальные веса вычисляются по формулам
w(x)=K \left( \frac{\rho(x,x_t)}{h(x_t)} \right ).
В этом случае отпадает необходимость задания параметра сглаживания f и его смысл эквивалентен выбору ширины окна h=h(x_t).

Робастные веса

Пусть

X^m\setminus\{x_t\} – обучающая выборка за исключением элемента x_i,
\hat{y_t}:=a \( x_t;X^m\setminus\{x_t\} \) – ответ алгоритма a, обученного на выборке X^m\setminus\{x_t\} при работе на объекте x_t.
\hat{\varepsilon_t}= \|\hat{y_t}-y_t \| – ошибка алгоритма на объекте x_i(ошибка скользящего контроля).

Пусть s - есть медиана величин\hat{\varepsilon_1},\ldots,\hat{\varepsilon_m}. тогда \delta_t = \bar{K}(\frac{\varepsilon_t}{6s}), где

\bar{K}(z)=\left{ (1-|z|^2)^2, \, \, \, |z|\le 1 \\ 0, \, \, \, |z|>1  \right.

Замечание

Возможны и другие варианты выбора весов \delta_t, например, занулить p штук, соответствующих наибольшим \hat{\varepsilon_t}. Это соответствует ядру
\bar{K}(z)=[z<\hat{\varepsilon}^{(m-p)}],

где \hat{\varepsilon}^{(m-p)} –- (m-p) - тый член вариационного ряда \hat{\varepsilon}^{(1)}<=\ldots<=\hat{\varepsilon}^{(m)}

В качестве весовой ядерной функции можно взять функцию Хубера (Huber, 1964) на которой основаны *[M-оценки]

K(z)=  \left{ z^2, \, \, \, |z| \le c\\ 2c|z|-c_2, \, \, \, |z|>c \right. Чтобы вычислить K(z) необходимо выбрать параметр устойчивости c. Одно популярное прикладное правило – c = 1,345 * s , где sробастная мера масштаба, такая как медианное абсолютное отклонение от медианы (MAD). Это популярное правило обеспечивает 95%-ую эффективность относительно гомоскедастичной нормальной модели в проблеме местоположения.

Алгоритм LOWESS

Вход

X^m - обучающая выборка;
w_t, \,\,\, t=1,\ldots,m весовые функции;

Выход

Коэффициенты \delta_t, \,\,\, t=1,\ldots,m

Алгоритм 1.1

1. Построить линейную регрессию во всех t=1,\ldots,t=m точках, используя весовые функции w_t, тем самым получим оценки для параметров модели \hat{\alpha_t}, \hat{\beta_t}.
А также приближения \hat{y_t}=\hat{\alpha_t}+\hat{\beta_t}x_t.
2. Инициализируем остатки \hat{\varepsilon_t}= \| \hat{y_t} - y_t \|. Вычислим робастные веса \delta_t =\bar{K}(\hat{\varepsilon_t})
3. повторять
4. Построить линейную регрессию во всех t=1,\ldots,t=m точках, используя весовые функции \delta_t w_t, тем самым получим оценки для параметров модели \hat{\alpha_t}, \hat{\beta_t}. А также приближения \hat{y_t}=\hat{\alpha_t}+\hat{\beta_t}x_t.
5. По новому набору значений \hat{\varepsilon_t}= \| \hat{y_t} - y_t \| вычислить новые значения коэффициентов \delta_t =\bar{K}(\hat{\varepsilon_t}).
6. пока веса \delta_t не стабилизируются

При использовании ядровых функций для оценки локальных весов объектов и робастных весов алгоритм модифицируется следующим образом:

Алгоритм 1.2

1. Инициализировать \delta_1:=\ldots=\delta_m:=1
2. повторять
3. Вычислить оценки скользящего контроля на каждом объекте
 \hat{y_t}:=a(x_t; X\setminus\{ x_t\}) = \frac{ \sum_{i=1, i\neq t }^{m} {y_i \delta_i K\left( \frac{\rho(x_i,x_t)}{h(x_t)}\right)} } {\sum_{i=1, i\neq t }^{m} {y_i K\left( \frac{\rho(x_i,x_t)}{h(x_t)}\right)} }
4. По набору значений \hat{\varepsilon_t}= \| \hat{y_t} - y_t \| вычислить новые значения коэффициентов \delta_t.
5. пока веса \delta_t не стабилизируются


Коэффициенты \delta_t, как и ошибки \varepsilon_t, зависят от функции a, которая, в свою очередь, зависит от \delta_t. На каждой итерации строится функция a, затем уточняются весовые множители \delta_t. Как правило, этот процесс сходится довольно быстро. Однако в практических реализациях имеет смысл вводить ограничение на количество итераций, как правило, это 2-3 итерации.

Примеры

Рис. 4. Пример для модельных данных
Рис. 4. Пример для модельных данных
На рисунке 4 представлен пример робастного локально-линейного сглаживания с помощь алгоритма LOWESS. С числом итераций цикла равным 2 и параметром сглаживания  f=0.5, то есть для приближения используется  k=\left\lfloor n/2 \right\rfloor ближайших точек выборки.

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

Литература

  1. A.I. McLeod Statistics 259b Robust Loess: S lowess. — 2004.
  2. Хардле В. Прикладная непараметрическая регрессия.. — Мир, 1993.
  3. Воронцов К.В. Лекции по алгоритмам восстановления регрессии. — 2007.
  4. John A Berger, Sampsa Hautaniemi, Anna-Kaarina Järvinen, Henrik Edgren, Sanjit K Mitra and Jaakko Astola Optimized LOWESS normalization parameter selection for DNA microarray data. — BMC Bioinformatics, 2004.
  5. Maronna, A., R. Martin, V. Yohai Robust Statistics: Theory and Methods.. — Wiley, 2006.
  6. Расин, Джеффри «Непараметрическая эконометрика: вводный курс». — Квантиль, №4, стр. 7–56., 2008.
  7. Huber, P.J. Robust estimation of a location parameter. Annals of Statistics 35. — 1964. — С. 73–101.

См. также

Личные инструменты