EM для смеси: мягкие назначения и граница гарантии

Если выводишь ответственности по Байесу, взвешенный M-шаг и неравенство монотонности, переходи к «Теперь сам».

Как думаешь, откуда появились два пика

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

Если бы метки были известны, мы разделили бы данные и оценили среднее и дисперсию каждого режима. Но меток нет, а жёстко назначать каждую точку ближайшему среднему слишком рано: области распределений могут перекрываться. EM сохраняет неопределённость назначения в виде весов.

На что это похоже

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

Это чередование гипотезы о происхождении и оценки параметров. Мягкая доля не равна произвольному коэффициенту похожести: она получается из конкретной смеси плотностей по формуле Байеса. Если меняется модель плотностей или ограничение параметров, шаги тоже нужно вывести заново.

Механизм: E-шаг и взвешенная цель

Обозначим доли α₁=γ, α₂=1−γ, где 0<γ<1, а нормальные плотности fⱼ имеют средние μⱼ и дисперсии vⱼ>0. При независимых наблюдениях лог-правдоподобие ℓ(θ)=Σᵢln(α₁f₁(xᵢ)+α₂f₂(xᵢ)). Параметр θ обозначает все доли, средние и дисперсии вместе.

E-шаг вычисляет ответственности rᵢⱼ=αⱼfⱼ(xᵢ)/(Σₖαₖfₖ(xᵢ)). Они неотрицательны и Σⱼrᵢⱼ=1. Формула - условная вероятность скрытого типа после наблюдения xᵢ: сначала берём совместную плотность значения и типа, затем нормируем по всем возможным типам.

На M-шаге r фиксированы. Максимизируется вспомогательная функция Q(θ|старое)=ΣᵢΣⱼrᵢⱼ[ln αⱼ+ln fⱼ(xᵢ)]. Это другая функция, а не прямая совместная максимизация исходной ℓ. В ней каждое наблюдение обучает каждую компоненту со своим весом. Положим Wⱼ=Σᵢrᵢⱼ, тогда ΣWⱼ=n.

Правило

При Wⱼ>0 M-шаг даёт αⱼ=Wⱼ/n, μⱼ=Σᵢrᵢⱼxᵢ/Wⱼ, vⱼ=Σᵢrᵢⱼ(xᵢ−μⱼ)²/Wⱼ. Дисперсия использует уже новое среднее и знаменатель Wⱼ, а не n−1. Для заданного ограничения vⱼ≥vmin>0 максимум этой же взвешенной цели имеет vⱼ=max(vmin,взвешенный средний квадрат). Произвольная поправка параметров не получает гарантию EM автоматически.

Средние и дисперсии выводятся тем же разложением квадратов, что для обычной MLE, только каждое слагаемое умножено на rᵢⱼ. Для долей сравни α с pⱼ=Wⱼ/n: разность цели равна nΣpⱼln(αⱼ/pⱼ) и не больше nΣpⱼ(αⱼ/pⱼ−1)=0 по неравенству ln z≤z−1. Поэтому доли тоже максимизируют Q при сумме 1.

Разбор: один M-шаг вручную

Наблюдения 0 и 2 получили ответственности первой компоненты 1/4 и 3/4. Её вес W₁=1, новая доля γ=1/2, среднее (0·1/4+2·3/4)/1=3/2. Дисперсия ((1/4)(3/2)²+(3/4)(1/2)²)/1=3/4. Вторая компонента имеет обратные веса, среднее 1/2 и ту же дисперсию 3/4.

Если применить невзвешенные формулы ко всем точкам, получатся одинаковые компоненты со средним 1 и дисперсией 1. Так потеряется именно то разделение, которое E-шаг передал M-шагу. Поэтому самостоятельная задача принимает дробные назначения в точном целочисленном масштабе и просит сохранить их веса.

import math

values = [-4, -3, -2, -1, 1, 2, 3, 4]


def e_step(data, params):
    gamma, m1, m2, v1, v2 = params
    responsibilities, total = [], 0.0
    for x in data:
        a = (math.log(gamma) - 0.5 * math.log(2 * math.pi * v1)
             - (x - m1) ** 2 / (2 * v1))
        b = (math.log(1 - gamma) - 0.5 * math.log(2 * math.pi * v2)
             - (x - m2) ** 2 / (2 * v2))
        top = max(a, b)
        left, right = math.exp(a - top), math.exp(b - top)
        responsibilities.append(left / (left + right))
        total += top + math.log(left + right)
    return responsibilities, total


def m_step(data, r, floor):
    weights = [sum(r), len(data) - sum(r)]
    if min(weights) <= 0:
        raise ValueError("компонента потерялась в машинном округлении")
    means, variances = [], []
    for j in range(2):
        w = r if j == 0 else [1 - p for p in r]
        mean = sum(p * x for p, x in zip(w, data)) / weights[j]
        variance = sum(p * (x - mean) ** 2
                       for p, x in zip(w, data)) / weights[j]
        means.append(mean)
        variances.append(max(floor, variance))
    return (weights[0] / len(data), *means, *variances)


params = (0.5, -1.0, 1.0, 4.0, 4.0)
for step in range(21):
    r, old = e_step(values, params)
    if step in [0, 1, 5, 20]:
        print("шаг", step, "ln L", f"{old:.6f}",
              "параметры", [round(x, 4) for x in params])
    new = m_step(values, r, 0.04)
    assert e_step(values, new)[1] >= old - 1e-10
    params = new
symmetric = (0.5, 0.0, 0.0, 7.5, 7.5)
r, old = e_step(values, symmetric)
print("симметрия сохраняется:", m_step(values, r, 0.04) == symmetric)
for variance in [0.01, 0.0001, 0.000001]:
    _, score = e_step([0, 1], (0.5, 0, 1, variance, 1))
    print("без ограничения v1 =", variance, "ln L =", f"{score:.6f}")

шаг 0 ln L -19.710503 параметры [0.5, -1.0, 1.0, 4.0, 4.0]
шаг 1 ln L -19.255298 параметры [0.5, -1.5302, 1.5302, 5.1584, 5.1584]
шаг 5 ln L -18.068686 параметры [0.5, -2.3244, 2.3244, 2.0972, 2.0972]
шаг 20 ln L -17.748674 параметры [0.5, -2.4887, 2.4887, 1.3066, 1.3066]
симметрия сохраняется: True
без ограничения v1 = 0.01 ln L = -0.862702
без ограничения v1 = 0.0001 ln L = 1.387046
без ограничения v1 = 1e-06 ln L = 3.684190

В распечатке параметры идут как γ,μ₁,μ₂,v₁,v₂. На первом шаге дисперсии выросли, а лог-правдоподобие улучшилось. Цель оценивает всю смесь, поэтому отдельный параметр не обязан монотонно расти или убывать. Программа выполняет фиксированное число шагов; строка 20 не является утверждением о доказанной сходимости. Числовая проверка допускает 1e−10 из-за округления.

Логарифмы плотностей a,b нормируются устойчиво: вычитаем max(a,b), экспоненты становятся не больше 1, затем возвращаем максимум к логарифму суммы. Это вычисление той же формулы, предотвращающее обнуление всех очень маленьких плотностей. При потере целой компоненты из-за машинного округления пример останавливается с сообщением, а не выдаёт деление на ноль за результат обучения.

Разбор: что именно доказывает монотонность

Для старых положительных ответственностей неравенство Йенсена для вогнутого логарифма даёт для каждого наблюдения ln Σⱼαⱼfⱼ(xᵢ) ≥ Σⱼrᵢⱼln(αⱼfⱼ(xᵢ)/rᵢⱼ). В правой части после суммирования по i получаются Q и постоянное по новым параметрам H=−Σᵢⱼrᵢⱼln rᵢⱼ.

На старых параметрах равенство достигается: отношения αⱼfⱼ(xᵢ)/rᵢⱼ одинаковы по j. M-шаг увеличивает Q, поэтому ℓ(новое)≥Q(новое|старое)+H≥Q(старое|старое)+H=ℓ(старое). Это нижняя касающаяся граница правдоподобия, а не утверждение, что каждый параметр по отдельности максимизирует исходное наблюдаемое ℓ при неизменных остальных.

Для нескольких компонент формулы те же с j=1,…,m. Для векторных наблюдений среднее считается покоординатно, а ковариация компоненты равна Σⱼ=Σᵢrᵢⱼ(xᵢ−μⱼ)(xᵢ−μⱼ)ᵀ/Wⱼ. Её положительную определённость проверяют отдельно: недостаточный разброс может дать вырожденную матрицу, к которой формула обычной нормальной плотности неприменима.

Что ломается без этого

Запускают обе компоненты с одинаковыми параметрами. Их ответственности тогда совпадают по всем наблюдениям; M-шаг сохраняет одинаковые средние и дисперсии. В симметричном примере программа показывает неподвижную точку с одинаковыми компонентами. Монотонность не заставляет алгоритм нарушить симметрию и найти более подходящее разделение. Разные старты могут давать разные результаты.

Ещё серьёзнее неограниченная смесь. Поставь среднее первой компоненты ровно в отдельное наблюдение и уменьшай её дисперсию к нулю, сохраняя положительную долю. Плотность в этой точке растёт как 1/√v₁. Вторая компонента даёт положительную плотность остальным наблюдениям. Поэтому всё правдоподобие может расти без предела. Последние строки программы демонстрируют такой путь на двух точках. Без ограничения дисперсий конечного глобального максимума может не быть.

Монотонность ℓ сама по себе не доказывает ни остановку за конечное число итераций, ни сходимость вектора параметров, ни обязательное попадание в локальный максимум. Условия компактности, регулярности и вид стационарных точек требуют отдельного анализа. Ограничение vmin делает показанный M-шаг корректной оптимизацией изменённой задачи, но не добавляет обещания глобальности.

Зачем это в работе

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

Теперь сам

Для трёх наблюдений ответственности первой компоненты равны 0,2;0,7;0,9. Какова новая γ? Что будет при одинаковых начальных компонентах? Обещает ли рост лог-правдоподобия конечный глобальный максимум для смеси без нижнего ограничения дисперсий?