Многомерная нормальная модель

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

До сих пор одному наблюдению соответствовало одно число $Y_i$. Во многих задачах компоненты отклика естественно моделировать совместно: несколько показателей одного объекта, измерения разных координат движения или связанные характеристики одного временного периода. Отдельные одномерные модели могут описать условные распределения каждой компоненты, но не задают их совместную зависимость.

Рассмотрим четырёхмерный нормальный отклик

$${\mathbf{Y}}_i=(Y_{i1},Y_{i2},Y_{i3},Y_{i4}{)}^{\mathsf{T}}\sim {\mathcal{N}}_4({\mathbf{\mu }}_i,\mathbf{\Sigma }),$$

где каждая компонента среднего линейно зависит от одного признака:

$${\mu }_{ij}={\beta }_{j,0}+{\beta }_{j,1}{x}_i,j=1,\dots ,4.$$

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

От коэффициентов к совместному распределению

Для каждой строки четыре блока среднего вычисляют ${\mathbf{\mu }}_i$. Четыре блока масштаба только со свободным членом задают положительные стандартные отклонения через логарифмическую link-функцию, а шесть блоков только со свободным членом задают частные корреляции на неограниченной шкале предиктора:

$${\sigma }_j=\exp ({\eta }_{\sigma ,j}),p_{jk}=\tanh ({\eta }_{p,jk}),k<j.$$

Путь вычисления имеет вид

$$\underset{18\text{коэффициентов}}{\underbrace{\mathbf{\beta }}}⟶\underset{4+4+6\text{предикторов}}{\underbrace{({\mathbf{\eta }}_{\mu ,i},{\mathbf{\eta }}_{\sigma },{\mathbf{\eta }}_p)}}⟶\underset{{\theta }_i}{\underbrace{({\mathbf{\mu }}_i,\mathbf{\sigma },\mathbf{p})}}⟶\mathbf{\Sigma }⟶{\mathcal{N}}_4({\mathbf{\mu }}_i,\mathbf{\Sigma }).$$

MvNormalMeanStdPartialCorrDefault<4> превращает частные корреляции в нижнетреугольный множитель корреляционной матрицы, а затем масштабирует его с помощью $\mathbf{\sigma }$. В результате

$$\mathbf{\Sigma }=\mathrm{diag}(\mathbf{\sigma })\mathbf{R}\mathrm{diag}(\mathbf{\sigma }),$$

где $\mathbf{R}$ является корректной положительно определённой корреляционной матрицей. Такая параметризация позволяет оптимизировать каждый предиктор на всей вещественной прямой и не проверять положительную определённость произвольной матрицы на каждой итерации.

Частная корреляция — не обычная парная корреляция

При $D=4$ строгая нижняя часть хранится по строкам:

$$(1,0),(2,0),(2,1),(3,0),(3,1),(3,2).$$

Первый параметр каждой строки совпадает с соответствующей обычной корреляцией, но последующие параметры имеют условный смысл. Например, $p_{21}$ описывает остаточную связь второй и первой компонент после учёта предшествующих компонент в выбранном порядке. Поэтому напечатанные partial_corr(2,1) и элемент $R_{21}$ корреляционной матрицы в общем случае различаются.

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

Воспроизводимая генерация данных

Синтетический пример задаёт четыре разные прямые средних, четыре маргинальных стандартных отклонения и шесть частных корреляций. Для каждого $x_i$ создаётся проверенный MvNormalMeanStdPartialCorrTheta, после чего try_fill_varying генерирует одну строку отклика для каждого набора параметров:

    const D: usize = 4;
    const N: usize = 400;
    const PARTIAL_CORRELATION_COUNT: usize = D * (D - 1) / 2;

    let x: [f64; N] = array::from_fn(|index| {
        let fraction = index as f64 / (N - 1) as f64;
        2.0_f64.mul_add(fraction, -1.0)
    });
    let true_mu_intercepts: [f64; D] = [1.0, -0.5, 0.25, 1.5];
    let true_mu_slopes: [f64; D] = [0.8, -0.6, 0.4, -0.3];
    let true_sigma: [f64; D] = [0.45, 0.70, 0.55, 0.80];
    // Packed row-major order: (1,0), (2,0), (2,1), (3,0), (3,1), (3,2).
    let true_partial_correlation_values = [0.45, -0.25, 0.30, 0.15, -0.35, 0.20];
    let true_partial_correlation =
        FixedPartialCorrelations::<D>::try_new(true_partial_correlation_values.to_vec())?;

    let family = MvNormalMeanStdPartialCorrDefault::<D>::new();
    let generating_theta = x
        .iter()
        .map(|&x_i| {
            let mu = array::from_fn(|component| {
                true_mu_slopes[component].mul_add(x_i, true_mu_intercepts[component])
            });
            MvNormalMeanStdPartialCorrTheta::try_new(mu, true_sigma, true_partial_correlation)
        })
        .collect::<Result<Vec<_>, _>>()?;
    let mut rng = StdRng::seed_from_u64(42);
    let mut y = vec![[0.0; D]; N];
    family.try_fill_varying(&mut rng, &generating_theta, &mut y)?;

Здесь используется то же ядро семейства, которое позднее вычисляет правдоподобие. Это исключает расхождение между формулами генератора и оцениваемой модели. TrySimulate возвращает SimulationError, если параметры недопустимы или численное преобразование породило недопустимое наблюдение; ошибки не кодируются значениями NaN.

Фиксированное зерно делает запуск воспроизводимым для закреплённых версий библиотек случайных чисел. При этом оценки не обязаны точно совпадать с истинными параметрами: они вычисляются по одной конечной случайной выборке.

Типизированные блоки многомерной модели

VectorParameterBlock<Mu, 4, ...> содержит четыре отдельных линейных предиктора с общим планом [1, x_i]. Аналогичный векторный блок для Sigma содержит четыре предиктора только со свободным членом. StrictLowerTriangularParameterBlock<PartialCorrelation, 4, ...> на уровне типа выражает отсутствие диагональных и верхнетреугольных коэффициентов зависимости.

    let mean_rows = x.iter().map(|&x_i| [1.0, x_i]).collect::<Vec<_>>();
    let mean_design = DenseDesign::from_rows(&mean_rows);
    let intercept = DenseDesign::intercept(N);
    let mu = VectorParameterBlock::<Mu, D, _, _>::new(
        array::from_fn(|_| LinearPredictorBlock::new(&mean_design)),
        NoPenalty,
        0,
    );
    let sigma = VectorParameterBlock::<Sigma, D, _, _>::new(
        array::from_fn(|_| LinearPredictorBlock::new(&intercept)),
        NoPenalty,
        mu.len(),
    );
    let partial_correlation =
        StrictLowerTriangularParameterBlock::<PartialCorrelation, D, _, _>::new(
            (0..PARTIAL_CORRELATION_COUNT)
                .map(|_| LinearPredictorBlock::new(&intercept))
                .collect(),
            NoPenalty,
            mu.len() + sigma.len(),
        );
    let blocks = ParameterBlocks::new((mu, sigma, partial_correlation));
    let mut model = Gamlss::try_new_with_observations(family, blocks, y.as_slice())?
        .with_objective_scale(ObjectiveScale::Mean);

Размеры блоков складываются следующим образом:

Конструктор ParameterBlocks::new назначает блокам последовательные участки глобального вектора. В примере смещения передаются и в отдельные конструкторы для наглядности, но итоговая совместимость всё равно проверяется при сборке модели.

В отличие от четырёх независимых нормальных моделей, один вызов правдоподобия использует всю строку [f64; 4] и общий ковариационный множитель. Поэтому производные по параметрам зависимости учитывают совместное положение всех четырёх компонент наблюдения.

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

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

На каждой итерации поиск шага уменьшает его, пока целевая функция не станет меньше. После успешного шага его размер осторожно увеличивается; цикл завершается, когда евклидова норма градиента становится меньше $10^{-6}$.

    let mut parameters = model.initial_parameters()?;
    let mut candidate = parameters.clone();
    let mut gradient = vec![0.0; model.nparams()];
    let mut step_size: f64 = 0.1;
    let mut iterations = 0;

    for iteration in 0..20_000 {
        model.gradient(&parameters, &mut gradient)?;
        let squared_gradient_norm = gradient.iter().map(|value| value * value).sum::<f64>();
        if squared_gradient_norm.sqrt() < 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.5);
                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={}",
                    squared_gradient_norm.sqrt()
                )
                .into());
            }
        }
        iterations = iteration + 1;
    }

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

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

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

    let diagnostics = model.training_diagnostics(&parameters)?;
    let coefficients = model.unpack_parameters(&parameters)?;
    let mu_coefficients = coefficients
        .blocks_of::<Mu>()
        .map(|block| [block.coefficients[0], block.coefficients[1]])
        .collect::<Vec<_>>();
    let fitted_theta = model.predict_theta_row(&parameters, 0)?;

    println!(
        "simple_multivariate_fit: n={N}, dimension={D}, parameters={}, iterations={iterations}, objective={:.6}, grad_norm={:.3e}",
        model.nparams(),
        diagnostics.objective,
        diagnostics.gradient_norm,
    );
    for component in 0..D {
        let observed = y
            .iter()
            .map(|observation| observation[component])
            .collect::<Vec<_>>();
        let pit = model.marginal_pit_values(&parameters, component, &observed)?;
        let (pit_min, pit_max) = finite_range(&pit);
        let pit_mean = finite_mean(&pit);
        println!(
            "component {component}: mu=[{:.4}, {:.4}] (true [{:.4}, {:.4}]), sigma={:.4} (true {:.4}), pit_mean={pit_mean:.4}, pit=[{pit_min:.4}, {pit_max:.4}]",
            mu_coefficients[component][0],
            mu_coefficients[component][1],
            true_mu_intercepts[component],
            true_mu_slopes[component],
            fitted_theta.sigma()[component],
            true_sigma[component],
        );
    }

    let mut packed_index = 0;
    for row in 1..D {
        for col in 0..row {
            let fitted = fitted_theta
                .partial_corr()
                .get(row, col)
                .expect("loop visits a strict-lower coordinate");
            println!(
                "partial_corr({row},{col})={fitted:.4} (true {:.4})",
                true_partial_correlation_values[packed_index],
            );
            packed_index += 1;
        }
    }

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

predict_theta_row возвращает полный набор параметров совместного распределения для одной строки. В примере достаточно первой строки, потому что sigma и частные корреляции имеют предикторы только со свободным членом; между строками меняется только mu. Если дать блокам масштаба или зависимости план с x, необходимо вычислять Theta для каждой новой строки.

Диагностика многомерной модели

Для каждой компоненты пример вычисляет маргинальные PIT

$$u_{ij}=\Phi \left(\frac{y_{ij}-{\mu }_{ij}}{{\sigma }_j}\right).$$

При корректных маргинальных распределениях значения $u_{ij}$ должны быть близки к равномерным. Напечатанные среднее и диапазон дают лишь быструю проверку; для содержательного анализа нужны гистограммы, эмпирические квантили и проверка зависимости PIT от признаков.

Маргинальные PIT не проверяют ковариационную структуру. Даже четыре идеально равномерных маргинальных PIT не отличат правильно оценённую зависимость от модели с независимыми компонентами. Для совместной диагностики семейство многомерных нормальных распределений поддерживает последовательное преобразование Розенблатта: каждая следующая компонента преобразуется через условную функцию распределения с учётом предыдущих. gamlss-diagnostics предоставляет rosenblatt_values, а интерпретация результата зависит от выбранного порядка компонент.

Ограничения и следующие расширения

В примере намеренно постоянны все параметры, кроме средних. Архитектура блоков позволяет независимо усложнять модель:

• дать отдельный план каждому Sigma и получить меняющуюся маргинальную неопределённость;
• моделировать выбранные частные корреляции через признаки;
• заменить линейные предикторы среднего сплайнами;
• перейти к MvStudentTCholesky, если совместные хвосты заметно тяжелее нормальных;
• использовать DynMvNormalCholesky, когда размерность известна только во время выполнения.

Число параметров зависимости растёт как $D(D-1)/2$. При большой размерности неструктурированная ковариационная модель быстро становится статистически и вычислительно дорогой; тогда нужны регуляризация, факторные или разреженные структуры и более сильные предположения о зависимости.