Двухкомпонентная гауссовская смесь для данных Old Faithful

Полный код примера находится в examples/faithful_mixture_fit.rs.

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

Рассмотрим классические данные Old Faithful. Одно наблюдение содержит длительность извержения в минутах и время ожидания следующего извержения в минутах:

$${\mathbf{Y}}_i=(Y_{i,\text{eruption}},Y_{i,\text{waiting}}{)}^{\mathsf{T}}.$$

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

$$p({\mathbf{y}}_i\mid \theta )=\sum _{c=0}^1{\pi }_c{\varphi }_2({\mathbf{y}}_i;{\mathbf{\mu }}_c,{\mathbf{\Sigma }}_c),{\pi }_c\ge 0,{\pi }_0+{\pi }_1=1.$$

Здесь $c$ обозначает скрытую компоненту смеси, а не координату многомерного отклика. У каждой из двух компонент есть двумерный вектор средних и собственная ковариационная матрица. Все предикторы в примере содержат только свободный член, поэтому параметры не зависят от номера наблюдения: это вероятностная кластеризация совместного распределения, а не регрессия по признакам.

Данные Old Faithful

Набор данных встроен в gamlss-datasets и загружается без операций ввода-вывода во время выполнения. data.x содержит номера наблюдений от 1 до 272, а data.y имеет тип &[[f64; 2]] и хранит пары [eruption, waiting]. Номера строк в этой модели не используются как предиктор.

    const D: usize = 2;
    const COMPONENTS: usize = 2;

    let data = gamlss_datasets::faithful();
    let y = data.y;
    let n = y.len();

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

Веса смеси и softmax с опорной компонентой

Веса должны принадлежать симплексу. Вместо независимой оптимизации двух ограниченных чисел модель использует один свободный логит $\alpha$ и фиксирует логит последней компоненты равным нулю:

$${\pi }_0=\frac{\exp (\alpha )}{\exp (\alpha )+1},{\pi }_1=\frac{1}{\exp (\alpha )+1}.$$

Такое softmax-представление с опорной компонентой автоматически даёт неотрицательные нормированные веса и устраняет лишнюю степень свободы. В коде ему соответствует SimplexLogitParameterBlock<MixtureWeight, 2, ...> с одним предиктором, содержащим только свободный член.

Для $C$ компонент нужны $C-1$ свободных логитов. Последняя компонента является опорной только в параметризации: это не делает её статистически особенной и не означает, что она является контрольной группой в данных.

Параметризация ковариационной матрицы

Каждая компонента задаётся семейством MvNormalCholeskyDefault<2>. Вместо непосредственной оптимизации трёх элементов симметричной ковариационной матрицы используется нижнетреугольный множитель Холецкого:

$${\mathbf{L}}_c=\left(\begin{array}{cc}\exp ({\eta }_{c,00})& 0\\ & \\ {\eta }_{c,10}& \exp ({\eta }_{c,11})\end{array}\right),{\mathbf{\Sigma }}_c={\mathbf{L}}_c{\mathbf{L}}_c^{\mathsf{T}}.$$

Логарифмические link-функции на диагонали гарантируют положительные диагональные элементы. Поэтому при конечных предикторах получается положительно определённая ковариационная матрица без отдельного учёта матричных ограничений при оптимизации. Её элементы равны

$${\mathbf{\Sigma }}_c=\left(\begin{array}{cc}{L}_{c,00}^2& L_{c,00}{L}_{c,10}\\ & \\ L_{c,00}{L}_{c,10}& L_{c,10}^2+L_{c,11}^2\end{array}\right).$$

Внедиагональный элемент $L_{c,10}$ сам по себе не является ни ковариацией, ни корреляцией. Для отчёта пример получает элементы ${\mathbf{\Sigma }}_c$ через theta.covariance и вычисляет

$${\rho }_c=\frac{{\Sigma }_{c,10}}{\sqrt{{\Sigma }_{c,00}{\Sigma }_{c,11}}}.$$

Вложенная структура параметров

Mixture<F, 2> повторяет семейство компоненты F два раза и добавляет перед ними симплекс весов. Для текущей модели структура имеет вид

$$(\text{weights},[({\mathbf{\mu }}_0,{\mathbf{L}}_0),({\mathbf{\mu }}_1,{\mathbf{L}}_1)]).$$

Эту же структуру выражают блоки: один SimplexLogitParameterBlock, а затем массив из двух пар VectorParameterBlock<Mu, 2, ...> и LowerTriangularParameterBlock<CholeskyScale, 2, ...>.

    let intercept = DenseDesign::intercept(n);
    let weights = SimplexLogitParameterBlock::<MixtureWeight, COMPONENTS, _, _>::new(
        vec![LinearPredictorBlock::new(&intercept)],
        NoPenalty,
        0,
    );
    let component_blocks = array::from_fn(|component| {
        let offset = weights.len() + component * 5;
        let mu = VectorParameterBlock::<Mu, D, _, _>::new(
            array::from_fn(|_| LinearPredictorBlock::new(&intercept)),
            NoPenalty,
            offset,
        );
        let cholesky = LowerTriangularParameterBlock::<CholeskyScale, D, _, _>::new(
            (0..3)
                .map(|_| LinearPredictorBlock::new(&intercept))
                .collect(),
            NoPenalty,
            offset + mu.len(),
        );
        (mu, cholesky)
    });
    let blocks = ParameterBlocks::new((weights, component_blocks));
    let family = Mixture::<_, COMPONENTS>::try_new(MvNormalCholeskyDefault::<D>::new())?;
    let mut model = Gamlss::try_new_with_observations(family, blocks, y)?
        .with_objective_scale(ObjectiveScale::Mean);

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

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

Правдоподобие и скрытая принадлежность

Отрицательное логарифмическое правдоподобие одного наблюдения равно

$${\ell }_i=-\log \left({\pi }_0{\varphi }_2({\mathbf{y}}_i;{\mathbf{\mu }}_0,{\mathbf{\Sigma }}_0)+{\pi }_1{\varphi }_2({\mathbf{y}}_i;{\mathbf{\mu }}_1,{\mathbf{\Sigma }}_1)\right).$$

Вычислять сумму плотностей напрямую численно ненадёжно, особенно в хвостах. Mixture объединяет логарифмические правдоподобия компонент с помощью приёма log-sum-exp, не вычисляя чрезвычайно малые плотности явно.

Апостериорная вероятность принадлежности наблюдения компоненте имеет вид

$$r_{ic}=P(Z_i=c\mid {\mathbf{y}}_i)=\frac{{\pi }_c{\varphi }_2({\mathbf{y}}_i;{\mathbf{\mu }}_c,{\mathbf{\Sigma }}_c)}{{\sum }_k{\pi }_k{\varphi }_2({\mathbf{y}}_i;{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)}.$$

Эти вероятности входят в аналитический градиент: производные по параметрам компоненты умножаются на её $r_{ic}$. На диаграмме рассеяния их можно показать непрерывной интенсивностью цвета либо назначать точке компоненту с максимальной вероятностью. Жёсткое назначение является производным представлением результата, а сама модель остаётся вероятностной.

Почему важна инициализация

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

Пример использует простую воспроизводимую эвристику: временно разделяет наблюдения по waiting < 70, а затем вычисляет внутри двух групп начальные средние, ковариационные матрицы и множители Холецкого. Начальный логит равен логарифму отношения размеров групп. Порог используется только для начального приближения; правдоподобие после этого оценивается на всех наблюдениях без известных меток компонент.

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

Оценивание параметров

Как и пример многомерной нормальной модели, эта глава сохраняет независимость библиотеки от оптимизатора и использует ручной градиентный спуск с поиском шага путём его последовательного сокращения. ObjectiveScale::Mean делит правдоподобие и его градиент на число наблюдений, чтобы характерный масштаб шага не зависел непосредственно от размера выборки.

    let starts = [
        empirical_component(y.iter().copied().filter(|row| row[1] < 70.0))?,
        empirical_component(y.iter().copied().filter(|row| row[1] >= 70.0))?,
    ];
    let mut parameters = initial_parameters(n, &starts);
    let mut candidate = parameters.clone();
    let mut gradient = vec![0.0; model.nparams()];
    let mut step_size: f64 = 0.05;
    let mut iterations = 0;

    for iteration in 0..20_000 {
        model.gradient(&parameters, &mut gradient)?;
        let gradient_norm = gradient
            .iter()
            .map(|value| value * value)
            .sum::<f64>()
            .sqrt();
        if gradient_norm < 1.0e-6 {
            iterations = iteration;
            break;
        }

        let objective = model.try_value(&parameters)?;
        loop {
            for ((next, parameter), gradient_value) in
                candidate.iter_mut().zip(&parameters).zip(&gradient)
            {
                *next = step_size.mul_add(-gradient_value, *parameter);
            }
            let candidate_objective = model.try_value(&candidate)?;
            if candidate_objective.is_finite() && candidate_objective < objective {
                parameters.copy_from_slice(&candidate);
                step_size = (step_size * 1.05).min(0.25);
                break;
            }
            step_size *= 0.5;
            if step_size < 1.0e-12 {
                return Err(format!(
                    "gradient descent line search failed at iteration {iteration}: objective={objective}, gradient_norm={gradient_norm}"
                )
                .into());
            }
        }
        iterations = iteration + 1;
    }

На каждой итерации пример принимает только шаг, уменьшающий целевую функцию, и завершает цикл при норме градиента меньше $10^{-6}$. Это прозрачный учебный оптимизатор, но не универсальный способ оценивания смесей: для прикладной модели важны несколько начальных приближений, явный статус сходимости, ограничения масштаба и сравнение найденных локальных оптимумов.

Запустить пример можно командой

cargo run --example faithful_mixture_fit --features multivariate

Интерпретация результата

Поскольку все блоки содержат только свободный член, predict_theta_row для любой строки возвращает один и тот же набор параметров на естественной шкале. Результат содержит нормированные веса и два MvNormalCholeskyTheta; для каждого из них доступны среднее, множитель Холецкого и элементы ковариационной матрицы.

    let diagnostics = model.training_diagnostics(&parameters)?;
    let fitted = model.predict_theta_row(&parameters, 0)?;

    println!(
        "faithful_mixture_fit: n={n}, components={COMPONENTS}, parameters={}, iterations={iterations}, objective={:.6}, grad_norm={:.3e}",
        model.nparams(),
        diagnostics.objective,
        diagnostics.gradient_norm,
    );
    for (component, (&weight, theta)) in
        fitted.weights().iter().zip(fitted.components()).enumerate()
    {
        let variance_eruption = theta.covariance(0, 0).expect("valid covariance index");
        let covariance = theta.covariance(1, 0).expect("valid covariance index");
        let variance_waiting = theta.covariance(1, 1).expect("valid covariance index");
        let correlation = covariance / (variance_eruption * variance_waiting).sqrt();
        println!(
            "component {component}: weight={weight:.4}, mean=[{:.4}, {:.4}], covariance=[[{variance_eruption:.4}, {covariance:.4}], [{covariance:.4}, {variance_waiting:.4}]], correlation={correlation:.4}",
            theta.mu()[0],
            theta.mu()[1],
        );
    }

Текущий детерминированный запуск даёт приблизительно следующие оценки:

Первая компонента соответствует режиму коротких извержений и короткого ожидания, вторая — режиму длинных извержений и долгого ожидания. Это содержательная интерпретация после оценивания, а не встроенное значение номера компоненты. При другой инициализации те же два режима могут получить противоположные индексы.

Для графика естественно показать исходную диаграмму рассеяния, оценённые средние и эллипсы постоянного расстояния Махаланобиса для каждой ${\mathbf{\Sigma }}_c$. Веса можно отразить в подписях, а вероятности принадлежности — цветом точек. Эллипсы описывают геометрию распределений компонент, но не являются жёсткими границами кластеров.

Идентифицируемость и вырождение

Перестановка компонент не меняет плотность смеси:

$${\pi }_0{\varphi }_0(\mathbf{y})+{\pi }_1{\varphi }_1(\mathbf{y})={\pi }_1{\varphi }_1(\mathbf{y})+{\pi }_0{\varphi }_0(\mathbf{y}).$$

Это переключение меток: параметры идентифицируются только с точностью до перестановки меток. Поэтому сравнивать компоненту 0 между двумя результатами оценивания можно лишь после осмысленного упорядочивания, например по среднему waiting.

У гауссовской смеси есть и более серьёзная особенность: правдоподобие не ограничено сверху, если одна ковариационная матрица стягивается к нулю вокруг отдельного наблюдения. На Old Faithful выбранное начальное приближение приводит к устойчивому содержательному решению, но это не устраняет проблему в общем случае. При прикладном оценивании полезны ограничения минимального масштаба, регуляризация ковариационных параметров, информативные априорные распределения или оптимизатор, специально учитывающий геометрию правдоподобия смеси.

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

Возможные расширения

Та же типизированная структура допускает несколько направлений развития:

• сделать веса смеси зависимыми от признаков и получить модель смеси экспертов;
• задать отдельные регрессионные предикторы для средних каждой компоненты;
• заменить двумерное нормальное распределение на MvStudentTCholesky, если внутри режимов нужны тяжёлые совместные хвосты;
• увеличить число компонент, изменив константный параметр типа C и добавив $C-1$ блоков логитов;
• добавить штрафы к параметрам компонент или коэффициентам механизма их выбора;
• генерировать новые наблюдения из оценённой смеси при включённой возможности rand.

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