Робастное оценивание и M-оценки

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

Точка разрыва и функция влияния

Точка разрыва (Breakdown Point, BP): Минимальная доля загрязнённых данных, при которой оценка ломается. X̄: BP = 1/n → 0. Медиана: BP = 0.5. S²: BP = 1/n. MAD = median|Xᵢ - median|: BP = 0.5.

Функция влияния (Hampel, 1974): IF(x; T, F) = lim_{ε→0} [T((1-ε)F + εδ_x) - T(F)]/ε. Для X̄: IF(x) = x - μ — неограничена. Для медианы: IF(x) = sign(x-μ)/(2f(μ)) — ограничена. Gross-error sensitivity: γ* = sup|IF(x)| < ∞ → робастная оценка.

M-оценки (Huber, 1964)

Определение: M-оценка θ̂ = argmin_θ Σᵢ ρ(Xᵢ - θ), или решение Σᵢ ψ(Xᵢ - θ) = 0, ψ = ρ'. Среднее: ρ(r) = r²/2. Медиана: ρ(r) = |r|. ОМП: ρ(r) = -log f(r).

M-оценка Хьюбера (c > 0): ρ_c(r) = r²/2 при |r| ≤ c; c|r| - c²/2 при |r| > c. ψ_c(r) = min(c, max(-c, r)). Среднее для малых остатков, медиана для больших. c = 1.345: 95% относительная эффективность при N(0,1).

Биквадратная оценка Тьюки: ψ(r) = r(1-(r/c)²)² при |r|≤c; 0 иначе. Полностью «отключает» экстремальные выбросы. c = 4.685: 95% эффективность. BP = 1/n (малая, но IF ограничена).

Алгоритм IRLS и асимптотика

IRLS (Iteratively Reweighted Least Squares): θ^{(t+1)} = Σwᵢ(θ^{(t)})Xᵢ / Σwᵢ(θ^{(t)}), wᵢ = ψ(rᵢ)/rᵢ. Каждая итерация — взвешенный МНК. Сходится к M-оценке. В линейной регрессии: β̂ = (XᵀWX)⁻¹XᵀWY.

Асимптотика M-оценок: √n(θ̂ - θ₀) →_d N(0, V), V = E[ψ²(X-θ₀)] / (E[ψ'(X-θ₀)])². Для оценки Хьюбера при N(0,1): V ≈ 1.08 (почти так же хорошо, как среднее!).

Оценка масштаба: Стандартизованный MAD: σ̂ = MAD/0.6745 — состоятельная оценка σ при нормальности, BP = 0.5.

Задание: (а) Данные: 1,2,3,4,5,50. X̄, медиана, усечённое среднее (20%), M-оценка Хьюбера (c=1.5). Как влияет выброс? (б) IRLS для линейной регрессии: данные с 2 выбросами. Сравните β̂_OLS и β̂_Huber — коэффициенты и их SE. (в) Симулируйте: n=100 из 0.9·N(0,1) + 0.1·N(10,1). Сравните дисперсию X̄ и медианы — какая оценка устойчивее?

Точка разрыва и устойчивость оценок

Точка разрыва (breakdown point) ε* — максимальная доля порчённых данных, при которой оценка остаётся ограниченной. X̄: ε* = 1/n ≈ 0 (один выброс может увести бесконечно). Медиана: ε* = 0.5 (наилучшая теоретически). Усечённое среднее (α-trimmed): ε* = α. Хьюбера: ε* зависит от c.

Влияние-функция (influence function): IF(x, T, F) = limₜ→₀ (T(F_t)−T(F))/t, где F_t = (1−t)F + tδ_x. Для среднего: IF(x,μ,F) = x−μ (неограниченная). Для медианы: IF(x,Med,F) = sgn(x−Med)/(2f(Med)) — ограниченная! Для M-оценки: IF(x,T,F) = ψ(x−T)/E[ψ'(X−T)].

Робастные оценки рассеяния и матрицы ковариации

MAD (Median Absolute Deviation): MAD = median|Xᵢ − median(X)|. Точка разрыва 50%. σ̂ = 1.4826·MAD (согласованность с нормальным).

Оценка MCD (Minimum Covariance Determinant): Для многомерных данных — найти подмножество h наблюдений с минимальным det(Cov). Точка разрыва ≈ (n−h)/n. Реализована в sklearn.covariance.MinCovDet. Применяется для выявления многомерных выбросов через расстояние Махаланобиса.

Тест на выбросы и методы обнаружения аномалий

Тест Грабса: H₀: нет выброса. Статистика G = max|Xᵢ−X̄|/s. Критическое значение зависит от n и α. Предполагает нормальность. Тест Тьюки: Выброс = точка за пределами Q₃+1.5IQR или Q₁−1.5IQR (ящик с усами). Непараметрический: не зависит от формы распределения. Для симметричных тяжёлых хвостов может ошибочно классифицировать.

Связь с теорией выживания

Функция выживания: S(t) = P(T > t). Функция риска: h(t) = f(t)/S(t) = −d/dt ln S(t). Для экспоненциального T: h(t) = λ (постоянный риск). Для Вейбулла: h(t) = λβ(λt)^{β−1} (убывает при β<1, возрастает при β>1). Оценка Каплана-Мейера: непараметрическая оценка S(t) при цензурировании. Ŝ(t) = Π_{tᵢ≤t}(1 − dᵢ/nᵢ), где dᵢ — число событий, nᵢ — число риска в момент tᵢ.

Копулы и зависимость хвостов

Копула C: [0,1]ⁿ → [0,1] — совместное распределение с равномерными маргиналами. Теорема Склара: F(x₁,...,xₙ) = C(F₁(x₁),...,Fₙ(xₙ)). Позволяет разделить моделирование маргиналей и зависимости. Гауссова копула: нет хвостовой зависимости (λ_L = λ_U = 0). Использовалась в CDO до 2008 — катастрофически занизила вероятность совместных дефолтов. Копула Клейтона: λ_L > 0 — хвостовая зависимость снизу. Применяется для кредитных рисков.

Многомерная статистика и анализ главных компонент

PCA: минимизировать ошибку реконструкции → собственные векторы матрицы ковариации. Свойство максимальной дисперсии: первые k главных компонент максимизируют дисперсию в k-мерном проекционном пространстве. Kernel PCA: замена Σ = XᵀX матрицей ядра K_ij = k(xᵢ,xⱼ) — нелинейная PCA. Sparse PCA (SPCA): добавить l₁-штраф к нагрузкам для интерпретируемости. Анализ независимых компонент (ICA): разложить смешанные сигналы на статистически независимые компоненты — «cocktail party problem».

Численный пример: метод главных компонент (PCA)

Задача: Матрица данных 4×2: точки (3,1),(1,3),(−3,−1),(−1,−3). Найти главные компоненты.

Шаг 1: x̄=(0,0). Матрица ковариаций: Σ = (1/4)·XᵀX. XᵀX=[(9+1+9+1, 3+3+3+3),(3+3+3+3, 1+9+1+9)]=[[20,12],[12,20]]. Σ=[[5,3],[3,5]].

Шаг 2: Собственные значения: det(Σ−λI)=(5−λ)²−9=0 → λ₁=8, λ₂=2.

Шаг 3: λ₁=8: вектор (1,1)/√2 (PC1 — диагональ). λ₂=2: вектор (1,−1)/√2 (PC2 — антидиагональ).

Шаг 4: PC1 объясняет 8/(8+2)=80% дисперсии. Проекция на PC1: z₁=(xᵢ₁+xᵢ₂)/√2 → {2√2, 2√2, −2√2, −2√2}≈{2.83, 2.83, −2.83, −2.83}. Сжатие с 2D до 1D с потерей лишь 20% дисперсии. ICA затем разделила бы проекцию на независимые источники сигнала.