Скользящее среднее по методу наименьших квадратов (МНК): локальная полиномиальная фильтрация данных без задержки и срезания пиков
Оглавление
- 1. Как работает локальная регрессия в скользящем окне
- 2. Почему расчет сводится к быстрой линейной свертке (метод Савицкого — Голея)
- 3. Таблица весовых коэффициентов и числовой пример расчета
- 4. Инструменты реализации: Python, Excel, R и Pine Script
- 5. Граничные эффекты: как обрабатывать края выборки
- 6. Вычисление производных без усиления шума
- 7. Граничные сценарии и частые сбои алгоритма
- 8. Как настроить параметры окна и полинома под реальный сигнал
- 9. Чек-лист проверки перед запуском в продакшен
- 10. План действий: от первых тестов к готовому конвейеру
- 11. Часто задаваемые вопросы (FAQ)
Простое скользящее среднее (SMA) заменяет значения в скользящем окне их средним арифметическим — горизонтальной константой. Из-за этого вершины пиков срезаются на 20–40%, а сам сигнал сдвигается вправо с фазовым запаздыванием на (N - 1) / 2 отсчетов.
Скользящее среднее по методу наименьших квадратов (МНК), известное в цифровой обработке сигналов как фильтр Савицкого — Голея, аппроксимирует точки каждого окна гибким полиномом (параболой или кубической кривой). Такой подход сохраняет форму, высоту и площадь пиков, обеспечивая нулевое запаздывание в центре окна при вычислительной скорости обычной линейной свертки.
| Метод сглаживания | Сохранение формы пиков | Фазовый сдвиг (запаздывание) | Вычислительная сложность | Основные сценарии применения |
|---|---|---|---|---|
| SMA (Простое скользящее) | Низкое: занижает амплитуду, размывает фронты | Высокий: сигнал сдвигается вправо на m отсчетов | O(1) при кольцевом буфере | Быстрая оценка общего тренда при отсутствии резких экстремумов |
| EMA (Экспоненциальное) | Среднее: искажает крутые и асимметричные всплески | Умеренный: снижен за счет весового акцента на свежие точки | Минимальная: один рекуррентный шаг | Сглаживание индикаторов в потоковом техническом анализе |
| МНК-фильтр (Савицкий — Голей) | Высокое: удерживает высоту, ширину и площадь импульса | Строго нулевой в центральном узле окна | Низкая: статическая линейная свертка с постоянным весовым ядром | Спектроскопия, хроматография, телеметрия датчиков, очистка признаков в ML |
1. Как работает локальная регрессия в скользящем окне
В основе традиционного сглаживания лежит скрытое допущение: внутри короткого интервала сигнал стабилен и может быть представлен горизонтальным отрезком. Скользящий метод наименьших квадратов снимает это ограничение, заменяя константу гибкой полиномиальной кривой.
| Параметр фрейма | Левое плечо окна | Центральный узел (отклик) | Правое плечо окна |
|---|---|---|---|
| Локальное время (t) | -m, …, -2, -1 | t = 0 | +1, +2, …, +m |
| Размер выборки | m отсчетов | 1 отсчет (y*0 = a0) | m отсчетов |
Процесс вычисления состоит из пяти шагов:
- Выделение симметричного окна. Для анализа берется отрезок нечетной длины:
N = 2m + 1
где m — полуширина окна. Нечетное число отсчетов гарантирует наличие единственной центральной точки x0, равноудаленной от левого (x-m … x-1) и правого (x1 … xm) краев выборки. - Центрирование локальной системы координат. Чтобы упростить расчеты, шкалу времени внутри окна сдвигают так, чтобы центральный узел получил координату t = 0. Левые точки нумеруются от -m до -1, правые — от +1 до +m.
- Построение локального полинома. В границах текущего фрейма строится регрессионная кривая степени k (чаще всего парабола k=2 или кубика k=3):
P(t) = ∑j=0k aj tj = a0 + a1t + a2t2 + … + aktk
Коэффициенты aj определяются классическим методом наименьших квадратов (OLS) через минимизацию суммы квадратов невязок между фактическими измерениями yi и значениями полинома P(ti):
S = ∑i=-mm (yi - P(ti))2 → min - Снятие отфильтрованного значения. Сглаженным значением для текущего положения окна становится точка полинома при t = 0:
y*0 = P(0) = a0 + a1·0 + a2·02 + … = a0
Все старшие степени переменной зануляются, поэтому отфильтрованная точка тождественна свободному коэффициенту a0. - Сдвиг окна. Окно смещается на один отсчет вправо, и процедура повторяется для следующей точки ряда.
2. Почему расчет сводится к быстрой линейной свертке (метод Савицкого — Голея)
Частая ошибка при первом знакомстве с методом — предположение, что на каждом шаге алгоритм заново формирует и обращает матрицы (XTX)-1XTy. В действительности всё считается намного быстрее.
Если шаг дискретизации ряда постоянен (Δt = const), локальные координаты отсчетов в любом окне строго зафиксированы: t ∈ {-m, …, -1, 0, 1, …, m}. Это означает, что матрица плана X и проекционная матрица (XTX)-1XT зависят исключительно от длины окна N и степени полинома k, оставаясь неизменными вдоль всего ряда данных.
В 1964 году Абрахам Савицкий и Марсель Голей доказали, что нахождение свободного члена a0 сводится к линейной комбинации отсчетов с постоянными целочисленными весами:
y*0 = (1 / H) ∑i=-mm ci · yi
где:
yi — исходные замеры внутри скользящего окна;
ci — целочисленные весовые коэффициенты ядра свертки;
H — нормировочный делитель, равный сумме всех весов ядра: H = ∑i=-mm ci.
Практическая закономерность: При расчете центрального узла окна весовые коэффициенты для полинома четной степени 2k полностью совпадают с коэффициентами для нечетной степени 2k+1. Параболическая регрессия (k=2) и кубическая регрессия (k=3) дают математически идентичный результат в центре окна. Разница между ними проявляется только при расчете производных или прогнозировании на краях диапазона.
3. Таблица весовых коэффициентов и числовой пример расчета
Коэффициенты симметричны относительно центрального элемента (c-i = ci). Это свойство обеспечивает строго нулевой фазовый сдвиг: пики и спады не смещаются во времени.
| Длина окна (N) | Степень полинома (k) | Нормировочный делитель (H) | Вектор весовых коэффициентов ядра (c-m … c0 … cm) | Готовая формула для MS Excel (для строки 5) |
|---|---|---|---|---|
| N=5 (m=2) | 2 или 3 | 35 | [-3, 12, 17, 12, -3] | =СУММПРОИЗВ(A3:A7; {-3;12;17;12;-3}) / 35 |
| N=7 (m=3) | 2 или 3 | 21 | [-2, 3, 6, 7, 6, 3, -2] | =СУММПРОИЗВ(A2:A8; {-2;3;6;7;6;3;-2}) / 21 |
| N=9 (m=4) | 2 или 3 | 231 | [-21, 14, 39, 54, 59, 54, 39, 14, -21] | =СУММПРОИЗВ(A1:A9; {-21;14;39;54;59;54;39;14;-21}) / 231 |
| N=11 (m=5) | 2 или 3 | 429 | [-36, 9, 44, 69, 84, 89, 84, 69, 44, 9, -36] | =СУММПРОИЗВ(A1:A11; {-36;9;44;69;84;89;84;69;44;9;-36}) / 429 |
Пошаговый расчет вручную
Разберем на реальных числах, почему МНК сохраняет высоту пика, в то время как простое усреднение его срезает.
Исходный фрагмент данных. Зафиксирован импульс на датчике давления из 5 точек (N=5, m=2):
y = [y-2 = 10, y-1 = 14, y0 = 25, y1 = 16, y2 = 11]
Фактический максимум находится в центре: y0 = 25.
Фильтрация простым средним (SMA):
SMA0 = (10 + 14 + 25 + 16 + 11) / 5 = 76 / 5 = 15.2
Результат: Вершина импульса упала с 25 до 15.2 (амплитуда срезана на 39.2%).
Фильтрация по МНК Савицкого — Голея (N=5, k=2, H=35):
Используем веса [-3, 12, 17, 12, -3]:
Sc = -3·10 + 12·14 + 17·25 + 12·16 + (-3)·11
Sc = -30 + 168 + 425 + 192 - 33 = 722
Нормирование суммы:
y*0 = Sc / H = 722 / 35 ≈ 20.63
Результат: Отфильтрованное значение составило 20.63. Отрицательные краевые множители (-3) компенсировали наклон боковых фронтов, предотвратив сплющивание вершины.
4. Инструменты реализации: Python, Excel, R и Pine Script
| Платформа / Библиотека | Основная область задач | Сильные стороны | Ограничения и риски |
|---|---|---|---|
| Python (scipy.signal) | Data Science, ML-конвейеры, инженерная обработка | Векторизованный C-код, поддержка граничных режимов, дифференцирование | Требует развертывания среды выполнения Python |
| MS Excel (формулы) | Экспресс-аудит небольших массивов, бизнес-отчеты | Наглядность, расчет без макросов и сторонних плагинов | Сложно масштабировать на таблицы объемом > 105 строк |
| R (signal, loess) | Биостатистика, эконометрика, научные публикации | Продвинутые тесты невязок, глубокая работа с локальной непараметрикой | Избыточен для простых конвейеров потоковой фильтрации |
| TradingView (Pine Script) | Сглаживание ценовых индикаторов | Нативная визуализация в торговом терминале | Риск эффекта заглядывания вперед (Look-ahead bias) на правом крае |
Реализация на Python (scipy.signal)
В экосистеме Python промышленным стандартом локальной МНК-фильтрации выступает функция savgol_filter:
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
# 1. Генерация тестового ряда с локальным экстремумом
np.random.seed(42)
t = np.linspace(0, 10, 200)
# Эталонный сигнал: гауссов импульс на линейном тренде
true_signal = np.exp(-((t - 5.0) ** 2) / 0.5) + 0.1 * t
# Наложение случайного шума
noisy_signal = true_signal + np.random.normal(0, 0.05, size=len(t))
# 2. Фильтрация методом Савицкого — Голея
# window_length: ширина окна (строго нечетное целое число)
# polyorder: степень полинома (2 — квадратичная, 3 — кубическая)
# mode: метод обработки границ ('interp' — локальная полиномиальная интерполяция)
filtered_signal = savgol_filter(
noisy_signal,
window_length=15,
polyorder=2,
mode='interp'
)5. Граничные эффекты: как обрабатывать края выборки
Для симметричного окна N = 2m + 1 крайние m точек слева и справа не имеют полного набора соседей. Например, при длине массива в 1000 элементов и окне N=11 (m=5) отсчеты с индексами от 0 до 4 и от 995 до 999 попадают в зону краевого дефицита.
| Стратегия обработки границы | Механика работы | Характер искажений | Сценарии применения |
|---|---|---|---|
| Усечение (drop) | Точки 0 … m-1 и L-m … L отбрасываются (заменяются на NaN) | Потеря 2m крайних отсчетов ряда | Оффлайн-анализ длинных исторических баз данных |
| Экстраполяция (mode='interp') | По крайнему полному окну строится единый полином, рассчитывающий крайние точки | Увеличение дисперсии ошибки на самых крайних 1–2 отсчетах | Научные замеры и выборки фиксированной длины |
| Зеркалирование (mode='mirror') | Сигнал зеркально отображается наружу относительно граничной координаты | Искусственно зануляет наклон производной на границе | Волновые колебания и циклические процессы |
| Асимметричные веса (Каузальный МНК) | Для каждой крайней точки строится отдельный вектор весов со смещением к центру | Возникновение фазового запаздывания на правом крае | Потоковые сигналы и системы реального времени |
Предупреждение о риске в потоковых данных (Look-ahead bias): В реальном времени симметричный фильтр нельзя применить к последней поступившей котировке без задержки на m тактов. Попытка использовать симметричную формулу на текущем баре означает использование еще не наступивших событий будущего. При тестировании торговой стратегии это создает иллюзию идеальной доходности, но на реальном счете приводит к систематическим убыткам из-за запаздывания одностороннего фильтра.
6. Вычисление производных без усиления шума
Анализ динамики процессов часто требует расчета скорости (первой производной) или ускорения (второй производной). Классическое численное дифференцирование через конечные разности приводит к взрывному росту высокочастотного шума:
(yi+1 - yi) / Δt
При делении случайной высокочастотной погрешности на малый шаг Δt амплитуда паразитных колебаний возрастает в десятки раз.
Так как скользящий МНК аппроксимирует точки окна непрерывной функцией P(t) = a0 + a1t + a2t2 + …, производные находятся прямым дифференцированием самого полинома в точке t = 0:
- Скорость (1-я производная): P'(0) = a1
- Ускорение (2-я производная): P''(0) = 2a2
На практике дифференцирование выполняется той же сверткой, но со специализированным дифференцирующим набором весовых коэффициентов ci(d) и обязательным масштабированием на шаг дискретизации в степени порядка производной Δtd:
y*0(d) = (d! / (Δtd · Hd)) ∑i=-mm ci(d) · yi
# Расчет сглаженной скорости (1-я производная, deriv=1)
velocity = savgol_filter(noisy_signal, window_length=11, polyorder=3, deriv=1, delta=0.05)
# Расчет сглаженного ускорения (2-я производная, deriv=2)
acceleration = savgol_filter(noisy_signal, window_length=11, polyorder=3, deriv=2, delta=0.05)7. Граничные сценарии и частые сбои алгоритма
Одиночные выбросы и импульсные помехи
Причина: Одиночный битый замер датчика (например, значение 1000 при фоновом уровне 15) квадратично штрафуется в методе наименьших квадратов и смещает всю локальную параболу на интервале окна N.
Решение: Предварительная очистка одномерным медианным фильтром малого порядка (scipy.ndimage.median_filter(signal, size=3)). Медиана устраняет спайки, а МНК-фильтр устраняет остаточный высокочастотный шум.
Осцилляции Рунге при завышении степени полинома
Причина: Задание степени k ≥ 4 при коротком окне (например, N=7, k=5) вызывает паразитные волны на краях интервала.
Решение: Для подавляющего большинства прикладных задач оптимальны порядки k=2 или k=3. Степени k ≥ 4 оправданы только в спектроскопии высокого разрешения при широких окнах (N ≥ 25).
Вырождение фильтра в интерполяцию (N ≤ k + 1)
Причина: Нулевое количество степеней свободы. Полином порядка k точно проходит через любые k+1 точек, повторяя весь исходный шум без сглаживания.
Решение: Всегда выдерживать достаточный запас по числу отсчетов: N ≥ k + 2.
8. Как настроить параметры окна и полинома под реальный сигнал
Выбор параметров N (ширины окна) и k (степени полинома) выполняется по следующему алгоритму:
- Оцените геометрию самого узкого пика. Найдите на графике информативный экстремум минимальной ширины и определите его полуширину в отсчетах (FWHM — Full Width at Half Maximum).
- Рассчитайте базовый размер окна:
Nбазовое ≈ 1.2 × FWHM
Условие: Значение N должно лежать в границах от 1.0 × FWHM до 1.5 × FWHM. Полученное число округлите вверх до ближайшего нечетного целого. - Выберите степень аппроксимации:
- ЕСЛИ ряд монотонный или содержит пологие участки тренда → выбирайте k=1 (локальная линейная регрессия).
- ЕСЛИ сигнал имеет выраженные вершины, провалы и точки перегиба → задавайте k=2 (парабола) или k=3 (кубический полином).
- Проведите анализ остатков (Residual Analysis):
Вычислите разность исходного и сглаженного рядов: ei = yi - y*i.
ЕСЛИ вектор e представляет собой случайный шум без выраженных трендов и автокорреляции → параметры подобраны верно.
ЕСЛИ в остатках прослеживаются контуры срезанных вершин полезного сигнала → окно N перераздуто, его следует уменьшить на 2–4 отсчета.
9. Чек-лист проверки перед запуском в продакшен
- [ ] Нечетный размер окна: Длина окна строго нечетная (N = 2m + 1 ∈ {5, 7, 9, 11, …}).
- [ ] Степени свободы: Соблюдено условие запаса точек N ≥ k + 2.
- [ ] Ограничение порядка полинома: Степень k ≤ 3, чтобы исключить краевые осцилляции Рунге.
- [ ] Защита от импульсных выбросов: Ряд проверен на единичные спайки; при их наличии настроен предварительный медианный фильтр.
- [ ] Корректный режим границ: Для оффлайн-рядов задан параметр mode='interp'; для real-time потоков заложена задержка выдачи на m тактов.
- [ ] Нормировка производных: При расчете производных сумма свертки поделена на Δtd.
- [ ] Равномерность сетки: Временной шаг Δt = const (при неравномерных интервалах свертка неприменима без интерполяции).
10. План действий: от первых тестов к готовому конвейеру
Шаг 1. Экспресс-тест в Excel (за 3 минуты)
- Вставьте столбец зашумленных данных в колонку A (начиная с ячейки A1).
- В ячейку B3 поместите формулу:
=СУММПРОИЗВ(A1:A5; {-3;12;17;12;-3}) / 35 - Растяните формулу вниз и постройте совместный график со столбцом простого усреднения
=СРЗНАЧ(A1:A5).
Шаг 2. Перенос в код на Python
Установите библиотеки:
pip install numpy scipy matplotlibПримените функцию savgol_filter к вашему датасету df['value'].values по примеру из Раздела 4.
Шаг 3. Тонкая калибровка
- Постройте разностный ряд e = y - y*.
- Скорректируйте размер окна N под ширину ваших экстремумов по правилу из Раздела 8.
