Введение

Эта книга — практическое введение в вероятностную регрессию и библиотеку gamlss. Она объясняет, как статистическая модель соотносится с типизированным API Rust, а исполняемые примеры в репозитории служат полным исходным кодом.

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

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

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

Вероятностная регрессия

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

Параметры условного распределения

Обычная линейная регрессия описала бы условное среднее наблюдений. Но неопределённость вокруг среднего тоже является частью предсказания: её можно увидеть как меняющийся разброс вокруг линии среднего.

Вероятностная регрессия позволяет модели описывать и среднее, и этот разброс.

Для семейства с K параметрами общая форма модели такова:

$$Y_i\mid x_i\sim D({\theta }_{i1},\dots ,{\theta }_{iK}),g_k({\theta }_{ik})={\eta }_{ik}.$$

У каждого параметра есть собственный предиктор. В линейном случае он имеет вид

$${\eta }_{ki}={\beta }_{k,0}+{\beta }_{k,1}{x}_{1i}+\dots +{\beta }_{k,p_k}{x}_{p_{k}i}.$$

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

Классы прикладных задач

Гетероскедастичность. Пусть средняя цена квартиры растёт с площадью, но разброс цен у больших квартир тоже выше. Тогда можно моделировать не только среднее $\mu$, но и масштаб $\sigma$ как функцию признаков.

Положительный ответ. Время ожидания, расходы и концентрации не могут быть отрицательными. Если нормальная модель иногда предсказывает невозможные отрицательные значения, уместнее выбрать $\mathrm{Gamma}$(гамма-распределение), $\mathrm{Weibull}$(распределение Вейбулла) или $\mathrm{LogNormal}$(логнормальное распределение).

Счётные данные и избыток нулей. Для числа обращений, покупок или дефектов естественны $\mathrm{Poisson}$(распределение Пуассона) и $\mathrm{NegBin}$(отрицательно-биномиальное распределение). Если нулей существенно больше, чем объясняет обычная счётная модель, может понадобиться модель с избыточными нулями.

Асимметрия и тяжёлые хвосты. Если данные заметно несимметричны или в них встречаются значения, слишком далёкие для нормального распределения, можно выбрать распределение с асимметрией или тяжёлыми хвостами: например, $\mathrm{SkewNormal}$(асимметричное нормальное) или распределение Стьюдента $t$.

Нормальная модель с линейным средним

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

Рассмотрим нормальную модель с линейным средним и постоянным стандартным отклонением:

$$Y_i\sim \mathcal{N}({\mu }_i,{\sigma }_i),$$

$${\mu }_i={\eta }_{\mu ,i}={\beta }_{\mu ,0}+{\beta }_{\mu ,1}{x}_i,{\sigma }_i=\exp ({\eta }_{\sigma ,i})=\exp ({\beta }_{\sigma ,0}).$$

В модели три глобальных коэффициента: свободный член и наклон для mu, а также свободный член на логарифмической шкале для sigma. Экспонента гарантирует положительность стандартного отклонения.

Вычисление параметров распределения

Для каждого наблюдения модель проходит один и тот же путь:

$$\underset{\text{3 глобальных коэффициента}}{\underbrace{\left(\begin{array}{c}{\beta }_{\mu ,0}\\ {\beta }_{\mu ,1}\\ {\beta }_{\sigma ,0}\end{array}\right)}}⟶\underset{\text{2 предиктора}}{\underbrace{\left(\begin{array}{c}{\eta }_{\mu ,i}\\ {\eta }_{\sigma ,i}\end{array}\right)}}⟶\underset{{\theta }_i}{\underbrace{\left(\begin{array}{c}{\mu }_i\\ {\sigma }_i\end{array}\right)}}⟶\mathcal{N}({\mu }_i,{\sigma }_i)⟶Y_i.$$

Первый переход выполняют блоки предикторов и матрицы плана, второй — link-функции, а семейство Normal вычисляет правдоподобие выбранного распределения.

Link-функции и блоки параметров

Предиктор eta живёт на удобной для оптимизации неограниченной шкале. Link-функция переводит его в допустимую область параметра распределения.

В нормальном примере используются:

$${\mu }_i={\eta }_{\mu ,i},{\sigma }_i=\exp ({\eta }_{\sigma ,i}).$$

Identity подходит для неограниченного параметра mu; Log подходит для положительного sigma. Ограничение положительности выражено и в типах: Normal требует, чтобы link-функция для Sigma реализовывала PositiveLink.

Структура блока параметра

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

В коде это разделение видно напрямую: для mu используется линейный план [1, x_i], а для постоянного sigma — только свободный член. Смещение mu.len() помещает коэффициенты sigma сразу после коэффициентов mu в глобальном векторе параметров.

    // 1. Assemble the typed model and its parameter blocks.
    let mu = ParameterBlock::<Mu, _, _>::linear(DenseDesign::from_rows(&x_design), NoPenalty, 0);
    let sigma =
        ParameterBlock::<Sigma, _, _>::linear(DenseDesign::intercept(n), NoPenalty, mu.len());
    let blocks = ParameterBlocks::new((mu, sigma));
    let mut model = Gamlss::try_new(Normal::<Identity, Log>::new(), blocks, &y)?;

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

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

$${\ell }_i=-\log p(y_i\mid {\theta }_i).$$

Оценивание минимизирует её сумму (и, при наличии, штрафы):

$$\mathcal{L}(\beta )=\sum _i-\log p(y_i\mid {\theta }_i(\beta ))+\mathrm{penalty}(\beta ).$$

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

В simple_fit.rs ради наглядности используется ручной градиентный спуск. Это учебный минимальный цикл, а не рекомендация для всех задач: для прикладных моделей необходимо выбрать оптимизатор, критерии остановки, масштабирование и обработку неудачных шагов.

    // 2. Fit the model with manual gradient-descent updates.
    let mut parameters = model.initial_parameters()?;
    let mut grad = vec![0.0; model.dim()];

    for _ in 0..10_000 {
        model.gradient(&parameters, &mut grad)?;
        for (parameter, grad_value) in parameters.iter_mut().zip(&grad) {
            *parameter -= 0.002 * grad_value;
        }
    }

Аналитический градиент

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

$$s_{\mu ,i}=\frac{\partial {\ell }_i}{\partial {\eta }_{\mu ,i}},s_{\sigma ,i}=\frac{\partial {\ell }_i}{\partial {\eta }_{\sigma ,i}}.$$

Затем блоки предикторов накапливают вклады наблюдений в глобальные коэффициенты:

$${\nabla }_{\beta }\mathcal{L}(\beta )=\sum _{i=1}^n\left(\begin{array}{c}{s}_{\mu ,i}\\ x_i{s}_{\mu ,i}\\ s_{\sigma ,i}\end{array}\right).$$

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

Оценённые параметры распределения и симуляция

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

В нормальной модели положения и масштаба это означает предсказывать mu_i и sigma_i; из них можно получить квантили, интервалы и вероятности событий.

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

    // 3. Predict distribution parameters and simulate fitted responses.
    let fitted_theta = model.predict_theta(&parameters)?;
    let fitted_rows = fitted_theta.len();

    #[cfg(feature = "rand")]
    let simulated = {
        use gamlss::core::TrySimulate;

        let mut simulation_rng = StdRng::seed_from_u64(43);
        let mut simulated = vec![f64::NAN; fitted_theta.len()];
        model
            .family()
            .try_fill_varying(&mut simulation_rng, &fitted_theta, &mut simulated)?;
        simulated
    };

Диагностика качества подгонки

Хорошая сходимость целевой функции ещё не означает, что распределение описывает данные адекватно. Для поддерживаемых семейств gamlss-diagnostics предоставляет значения функции распределения (CDF) и PIT, нормированные квантильные остатки и CRPS.

В simple_fit.rs доступны методы pit_values и quantile_residuals. Для корректной модели PIT-значения в совокупности должны быть близки к равномерному распределению, а нормализованные квантильные остатки — к стандартному нормальному.

    // 4. Inspect diagnostics for the fitted model.
    let diagnostics = model.training_diagnostics(&parameters)?;
    let pit = model.pit_values(&parameters)?;
    let residuals = model.quantile_residuals(&parameters)?;

Ограничения базовой модели

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

В следующей модели семейство и предиктор mu останутся прежними, а план для Sigma изменится с [1] на [1, x_i]. Такое минимальное изменение покажет ситуацию, в которой RMSE почти не улучшается, а вероятностный прогноз становится заметно лучше.

Нормальная модель с переменным масштабом

Это следующая модель после simple_fit: она меняет только один компонент уже знакомой нормальной модели и впервые показывает вероятностную регрессию в действии.

Постановка задачи

Сгенерируем достаточно много наблюдений на равномерной сетке $x\in [-2,2]$ по модели

$$\mu (x)=1+0.8x,\sigma (x)=\exp (-0.7+0.45x),Y\mid x\sim \mathcal{N}(\mu (x),\sigma (x)).$$

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

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

Сравниваемые спецификации модели

Базовая модель повторяет предыдущую главу:

$$M_0:{\mu }_i={\beta }_0+{\beta }_1{x}_i,\log {\sigma }_i={\gamma }_0.$$

Вероятностная модель добавляет x только в предиктор масштаба:

$$M_1:{\mu }_i={\beta }_0+{\beta }_1{x}_i,\log {\sigma }_i={\gamma }_0+{\gamma }_1{x}_i.$$

У $M_0$ три коэффициента, у $M_1$ четыре. Семейство остаётся Normal::<Identity, Log>, link-функции остаются прежними, а блок mu вообще не меняется. В коде Rust существенное изменение сводится к замене плана только со свободным членом для Sigma на план с колонками [1, x_i]:

let sigma_rows = x.map(|x_i| [1.0, x_i]);
let sigma = ParameterBlock::<Sigma, _, _>::linear(
    DenseDesign::from_rows(&sigma_rows),
    NoPenalty,
    mu.len(),
);

Именно это небольшое изменение составляет смысл главы: в GAMLSS новый регрессионный вопрос обычно выражается новым блоком предиктора для конкретного параметра, а не заменой всей модели.

Критерии сравнения

Сравнение должно быть проведено на одних и тех же тестовых строках и включать не только точечную метрику.

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

Интерпретация коэффициента предиктора масштаба

Так как используется логарифмическая link-функция,

$$\frac{\sigma (x+1)}{\sigma (x)}=\exp ({\gamma }_1).$$

Следовательно, gamma_1 задаёт мультипликативное изменение условного стандартного отклонения при увеличении x на единицу. Это не эффект x на средний отклик и не добавка к дисперсии. Глава должна показать интерпретацию и на шкале предиктора, и на естественной шкале sigma.

Границы рассматриваемой модели

В этой главе не нужны новое семейство, сплайны, штрафы, категориальные признаки, преобразование отклика или сложный API оптимизатора. Надёжный оптимизатор можно подключить через уже готовый адаптер, но его устройство следует вынести в отдельную техническую вставку. Одновременное добавление нескольких возможностей уничтожит контролируемое сравнение $M_0$ и $M_1$.

Структура воспроизводимого примера

Полный examples/heteroscedastic_normal.rs должен выполнять одну воспроизводимую последовательность действий:

1. Сгенерировать данные и выполнить фиксированное чередующееся разбиение.
2. Собрать $M_0$ и $M_1$ из одинакового плана mu.
3. Обучить обе модели с одинаковыми критериями остановки.
4. Проверить конечность целевой функции, норму градиента и статус оптимизатора.
5. Получить параметры, квантили и интервалы на тестовой части.
6. Напечатать одну компактную таблицу с RMSE, NLL, CRPS и покрытием.
7. Сохранить числовые ряды, из которых воспроизводится иллюстрация для книги.

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

После этой главы переход к нелинейным mu(x) и sigma(x) требует заменить линейные матрицы плана на сплайн-предикторы, но не требует заново объяснять, зачем у масштаба есть собственная регрессия.

Сплайн-предикторы положения, масштаба и асимметрии

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

Предыдущие модели поочерёдно усложняли нормальное условное распределение: сначала признаки объясняли только линейное среднее, затем собственный предиктор появился у масштаба. Теперь добавим нелинейные предикторы, явную асимметрию распределения и прикладной оптимизатор с поиском шага. Каждому параметру семейства по-прежнему соответствует отдельный типизированный блок.

Данные a1

Встроенный набор данных содержит 2000 пар $(x_i,y_i)$ с $x\in [0,1.2]$. Он загружается из gamlss-datasets без операций ввода-вывода во время выполнения и дополнительных выделений памяти:

    let data = gamlss_datasets::a1();
    let (x, y) = (data.x, data.y);
    let n = y.len();

Зависимость среднего от x явно нелинейна. Ширина облака также меняется по диапазону признака, а форма отклонений не везде симметрична. На рисунке показаны результат оценивания и отдельная кривая параметра асимметрии:

Отклик принимает отрицательные значения, поэтому гамма- и логнормальное распределения на исходной шкале не подходят. Асимметричное нормальное распределение определено на всей вещественной прямой и позволяет моделировать наблюдаемую асимметрию без преобразования y.

Почему параметризация через среднее и стандартное отклонение

В примере используется SkewNormalMeanSd, параметры которой на естественной шкале имеют прямой смысл:

$$Y_i\mid x_i\sim {\mathrm{SN}}_{\text{mean/sd}}(m_i,s_i,{\nu }_i),\mathbb{E}[Y_i\mid x_i]=m_i,\mathrm{SD}(Y_i\mid x_i)=s_i.$$

Параметр ${\nu }_i$ управляет направлением и силой асимметрии; при ${\nu }_i=0$ семейство совпадает с нормальным распределением. Параметризация через среднее и стандартное отклонение выбрана потому, что первый сплайн остаётся условным средним даже при меняющемся $\nu (x)$.

Link-функции в модели имеют вид

$$m_i={\eta }_{m,i},s_i=\exp \left(\mathrm{clamp}({\eta }_{s,i},-8,2)\right),{\nu }_i={\eta }_{\nu ,i}.$$

ClampedLog<-8, 2> гарантирует положительное стандартное отклонение и ограничивает экстремальные значения масштаба во время оптимизации. В найденном решении ограничения не активны.

Три сплайн-предиктора

Для каждого параметра строится собственный базис кубических B-сплайнов:

$${\eta }_r(x)=\sum _{j=0}^{K_r-1}{\beta }_{r,j}{B}_{r,j}(x),r\in \{m,s,\nu \}.$$

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

Всего оптимизатор видит $20+14+12=46$ коэффициентов.

    // Each spline contains its own constant direction because the B-spline
    // rows sum to one. Second-difference penalties suppress local wiggles while
    // leaving constant and linear trends unpenalized.
    let mean_design = OpenUniformSplineDesign::from_data(x, MEAN_BASIS, SplineOrder::Cubic)?;
    let sigma_design = OpenUniformSplineDesign::from_data(x, SIGMA_BASIS, SplineOrder::Cubic)?;
    let nu_design = OpenUniformSplineDesign::from_data(x, NU_BASIS, SplineOrder::Cubic)?;

    let mean = ParameterBlock::<Mean, _, _>::new(
        mean_design,
        PreparedDifferencePenalty::try_new(MEAN_SMOOTHING, 2)?,
        0,
    );
    let sigma_offset = mean.len();
    let sigma = ParameterBlock::<Sigma, _, _>::new(
        sigma_design,
        PreparedDifferencePenalty::try_new(SIGMA_SMOOTHING, 2)?,
        sigma_offset,
    );
    let nu_offset = sigma_offset + sigma.len();
    let nu = ParameterBlock::<Nu, _, _>::new(
        nu_design,
        PreparedDifferencePenalty::try_new(NU_SMOOTHING, 2)?,
        nu_offset,
    );

    // The mean/SD parameterization keeps the first two fitted curves directly
    // interpretable even when the skewness parameter changes with x.
    let model = Gamlss::try_new(
        SkewNormalMeanSd::<Identity, ClampedLog<-8, 2>, Identity>::new(),
        ParameterBlocks::new((mean, sigma, nu)),
        y,
    )?
    .with_objective_scale(ObjectiveScale::Mean);

Смещения вычисляются из len() предыдущих блоков. В результате структура плоского вектора коэффициентов остаётся явной, а соответствие (Mean, Sigma, Nu) требованиям семейства проверяется при компиляции и при построении модели.

Штрафы P-сплайнов

Большой базис сам по себе только даёт модели возможность изгибаться. Гладкость задаёт PreparedDifferencePenalty второго порядка. Для одного блока с $K$ коэффициентами библиотека добавляет

$$J_2(\mathbf{\beta })=\frac{\lambda }{K-2}\sum _{j=0}^{K-3}{\left({\beta }_{j+2}-2{\beta }_{j+1}+{\beta }_j\right)}^2.$$

Комбинацию базиса B-сплайнов и штрафа на разности обычно называют P-сплайном. Штраф подавляет резкие изменения наклона соседних коэффициентов, но не штрафует их постоянный и линейный профили.

В этом примере числа базисных функций и значения $\lambda$ зафиксированы вручную.

Модель использует ObjectiveScale::Mean, поэтому оптимизируется

$$\mathcal{L}(\mathbf{\beta })=\frac{1}{n}\sum _{i=1}^n-\log p(y_i\mid m_i,s_i,{\nu }_i)+J_m+J_s+J_{\nu }.$$

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

Инициализация и два старта

Строки базиса B-сплайнов образуют разбиение единицы. Поэтому одинаковые коэффициенты одного блока задают константу для всех наблюдений. Начальные коэффициенты сплайна среднего равны общему выборочному среднему, а сплайна логарифма масштаба — логарифму общего стандартного отклонения. Блок асимметрии запускается из двух констант: $\nu =-5$ и $\nu =5$.

Ненулевой старт здесь принципиален. В параметризации через среднее и стандартное отклонение $\nu =0$ является стационарной точкой симметричной нормальной подмодели. Градиентный оптимизатор, запущенный ровно из нуля, может остаться в ней даже тогда, когда решение с асимметрией имеет меньшее отрицательное логарифмическое правдоподобие. Два знака не гарантируют глобального минимума, но устраняют наиболее очевидную зависимость от этой симметричной точки; пример сохраняет решение с меньшей целевой функцией со штрафом.

    // A constant start on the natural response scale is more useful than an
    // all-zero spline start. Equal B-spline coefficients produce a constant.
    let (sample_mean, sample_sd) = mean_sd(y);
    let mut constant_start = model.initial_parameters()?;
    constant_start[..sigma_offset].fill(sample_mean);
    constant_start[sigma_offset..nu_offset].fill(sample_sd.ln());

    // The adapter retains the reusable GAMLSS workspace behind RefCell because
    // argmin asks for `&self`, whereas Objective deliberately uses `&mut self`.
    // nu=0 is a stationary symmetric submodel in the mean/SD
    // parameterization. We therefore fit both non-zero signs and keep the
    // solution with the lower penalized objective.
    let nu_range = nu_offset..model.nparams();
    let candidate_fits = NU_STARTS
        .into_iter()
        .map(|nu_start| {
            fit_once(
                model.clone().into_workspace_objective(),
                &constant_start,
                nu_range.clone(),
                nu_start,
            )
        })
        .collect::<Result<Vec<_>, _>>()?;

Значения $\pm 5$ используются только как начальные предикторы асимметрии; L-BFGS свободно изменяет все 12 коэффициентов Nu.

Подключение argmin

gamlss-core намеренно не зависит от конкретного оптимизатора. Его Objective предоставляет значение и аналитический градиент над обычным &[f64], а тонкий адаптер в примере реализует типажи CostFunction и Gradient из argmin.

/// Single-threaded adapter from the mutable, buffer-reusing GAMLSS objective to
/// the shared-reference callbacks expected by argmin.
#[derive(Debug)]
struct ArgminObjective<O> {
    objective: RefCell<O>,
}

impl<O> ArgminObjective<O> {
    const fn new(objective: O) -> Self {
        Self {
            objective: RefCell::new(objective),
        }
    }

    fn dim(&self) -> usize
    where
        O: Objective,
    {
        self.objective.borrow().dim()
    }
}

impl<O> CostFunction for ArgminObjective<O>
where
    O: Objective,
    O::Error: std::error::Error + Send + Sync + 'static,
{
    type Param = Vec<f64>;
    type Output = f64;

    fn cost(&self, param: &Self::Param) -> Result<Self::Output, Error> {
        self.objective.borrow_mut().value(param).map_err(Error::new)
    }
}

impl<O> Gradient for ArgminObjective<O>
where
    O: Objective,
    O::Error: std::error::Error + Send + Sync + 'static,
{
    type Param = Vec<f64>;
    type Gradient = Vec<f64>;

    fn gradient(&self, param: &Self::Param) -> Result<Self::Gradient, Error> {
        let mut gradient = vec![0.0; self.dim()];
        self.objective
            .borrow_mut()
            .gradient(param, &mut gradient)
            .map_err(Error::new)?;
        Ok(gradient)
    }
}
fn main() -> Result<(), Box<dyn std::error::Error>> {
    use gamlss::core::{
        ClampedLog, Gamlss, HasQuantile, Identity, Mean, Nu, ObjectiveScale, ParameterBlock,
        ParameterBlocks, Sigma,
    };
    use gamlss::diagnostics::CdfDiagnosticsExt;
    use gamlss::family::SkewNormalMeanSd;
    use gamlss::spline::{OpenUniformSplineDesign, PreparedDifferencePenalty, SplineOrder};

    let output_mode = OutputMode::from_args()?;
    let data = gamlss_datasets::a1();
    let (x, y) = (data.x, data.y);
    let n = y.len();

    // Each spline contains its own constant direction because the B-spline
    // rows sum to one. Second-difference penalties suppress local wiggles while
    // leaving constant and linear trends unpenalized.
    let mean_design = OpenUniformSplineDesign::from_data(x, MEAN_BASIS, SplineOrder::Cubic)?;
    let sigma_design = OpenUniformSplineDesign::from_data(x, SIGMA_BASIS, SplineOrder::Cubic)?;
    let nu_design = OpenUniformSplineDesign::from_data(x, NU_BASIS, SplineOrder::Cubic)?;

    let mean = ParameterBlock::<Mean, _, _>::new(
        mean_design,
        PreparedDifferencePenalty::try_new(MEAN_SMOOTHING, 2)?,
        0,
    );
    let sigma_offset = mean.len();
    let sigma = ParameterBlock::<Sigma, _, _>::new(
        sigma_design,
        PreparedDifferencePenalty::try_new(SIGMA_SMOOTHING, 2)?,
        sigma_offset,
    );
    let nu_offset = sigma_offset + sigma.len();
    let nu = ParameterBlock::<Nu, _, _>::new(
        nu_design,
        PreparedDifferencePenalty::try_new(NU_SMOOTHING, 2)?,
        nu_offset,
    );

    // The mean/SD parameterization keeps the first two fitted curves directly
    // interpretable even when the skewness parameter changes with x.
    let model = Gamlss::try_new(
        SkewNormalMeanSd::<Identity, ClampedLog<-8, 2>, Identity>::new(),
        ParameterBlocks::new((mean, sigma, nu)),
        y,
    )?
    .with_objective_scale(ObjectiveScale::Mean);

    // A constant start on the natural response scale is more useful than an
    // all-zero spline start. Equal B-spline coefficients produce a constant.
    let (sample_mean, sample_sd) = mean_sd(y);
    let mut constant_start = model.initial_parameters()?;
    constant_start[..sigma_offset].fill(sample_mean);
    constant_start[sigma_offset..nu_offset].fill(sample_sd.ln());

    // The adapter retains the reusable GAMLSS workspace behind RefCell because
    // argmin asks for `&self`, whereas Objective deliberately uses `&mut self`.
    // nu=0 is a stationary symmetric submodel in the mean/SD
    // parameterization. We therefore fit both non-zero signs and keep the
    // solution with the lower penalized objective.
    let nu_range = nu_offset..model.nparams();
    let candidate_fits = NU_STARTS
        .into_iter()
        .map(|nu_start| {
            fit_once(
                model.clone().into_workspace_objective(),
                &constant_start,
                nu_range.clone(),
                nu_start,
            )
        })
        .collect::<Result<Vec<_>, _>>()?;

    let fit = candidate_fits
        .into_iter()
        .min_by(|left, right| left.objective.total_cmp(&right.objective))
        .expect("NU_STARTS is non-empty");
    let parameters = fit.parameters;

    let fitted = model.predict_theta(&parameters)?;
    let intervals_90 = fitted
        .iter()
        .map(|theta| {
            (
                model.family().quantile(0.05, theta),
                model.family().quantile(0.95, theta),
            )
        })
        .collect::<Vec<_>>();
    if output_mode == OutputMode::Csv {
        print_prediction_csv(x, y, &fitted, &intervals_90);
        return Ok(());
    }

    let diagnostics = model.training_diagnostics(&parameters)?;
    let pit = model.pit_values(&parameters)?;
    let residuals = model.quantile_residuals(&parameters)?;
    let coverage_90 = intervals_90
        .iter()
        .zip(y)
        .filter(|(interval, observation)| (interval.0..=interval.1).contains(observation))
        .count() as f64
        / n as f64;

    let (mean_min, mean_max) = range(fitted.iter().map(|theta| theta.mean));
    let (sigma_min, sigma_max) = range(fitted.iter().map(|theta| theta.sigma));
    let (nu_min, nu_max) = range(fitted.iter().map(|theta| theta.nu));
    let (pit_mean, pit_sd) = mean_sd(&pit);
    let (residual_mean, residual_sd) = mean_sd(&residuals);

    println!(
        "a1 skew-normal spline fit: n={n}, parameters={}, starts={}, selected_nu_start={:.1}",
        parameters.len(),
        NU_STARTS.len(),
        fit.nu_start,
    );
    println!(
        "iterations={}, termination={}",
        fit.iterations, fit.termination,
    );
    println!(
        "objective={:.6}, mean_nll={:.6}, penalty={:.6}, grad_norm={:.3e}",
        diagnostics.objective,
        diagnostics.train_nll,
        diagnostics.penalty,
        diagnostics.gradient_norm,
    );
    println!(
        "fitted mean=[{mean_min:.4}, {mean_max:.4}], sigma=[{sigma_min:.4}, {sigma_max:.4}], nu=[{nu_min:.4}, {nu_max:.4}]",
    );
    println!(
        "PIT mean={pit_mean:.4}, PIT sd={pit_sd:.4}, qres mean={residual_mean:.4}, qres sd={residual_sd:.4}, 90% coverage={coverage_90:.3}",
    );
    Ok(())
}

fn mean_sd(values: &[f64]) -> (f64, f64) {
    debug_assert!(!values.is_empty());
    let count = values.len() as f64;
    let mean = values.iter().sum::<f64>() / count;
    let variance = values
        .iter()
        .map(|value| {
            let centered = value - mean;
            centered * centered
        })
        .sum::<f64>()
        / count;
    (mean, variance.sqrt())
}

fn range(mut values: impl Iterator<Item = f64>) -> (f64, f64) {
    let first = values.next().expect("fitted values are non-empty");
    values.fold((first, first), |(min, max), value| {
        (min.min(value), max.max(value))
    })
}

fn fit_once<O>(
    objective: O,
    constant_start: &[f64],
    nu_range: Range<usize>,
    nu_start: f64,
) -> Result<Fit, Error>
where
    O: Objective,
    O::Error: std::error::Error + Send + Sync + 'static,
{
    let mut initial_parameters = constant_start.to_vec();
    initial_parameters[nu_range].fill(nu_start);

    let problem = ArgminObjective::new(objective);
    let line_search = MoreThuenteLineSearch::new().with_c(1.0e-4, 0.9)?;
    let solver = LBFGS::new(line_search, LBFGS_MEMORY)
        .with_tolerance_grad(1.0e-5)?
        .with_tolerance_cost(1.0e-8)?;
    let result = Executor::new(problem, solver)
        .configure(|state| state.param(initial_parameters).max_iters(MAX_ITERATIONS))
        .run()?;

    let parameters = result
        .state
        .get_best_param()
        .or_else(|| result.state.get_param())
        .cloned()
        .ok_or_else(|| Error::msg("argmin finished without a parameter vector"))?;

    Ok(Fit {
        parameters,
        objective: result.state.get_best_cost(),
        iterations: result.state.get_iter(),
        termination: result
            .state
            .get_termination_reason()
            .map_or_else(|| "unknown".to_owned(), ToString::to_string),
        nu_start,
    })
}

fn print_prediction_csv(
    x: &[f64],
    y: &[f64],
    fitted: &[gamlss::family::SkewNormalMeanSdTheta],
    intervals_90: &[(f64, f64)],
) {
    debug_assert_eq!(x.len(), y.len());
    debug_assert_eq!(x.len(), fitted.len());
    debug_assert_eq!(x.len(), intervals_90.len());
    println!("x,y,mean,sigma,nu,q05,q95");
    for (((x, observation), theta), &(q05, q95)) in x.iter().zip(y).zip(fitted).zip(intervals_90) {
        println!(
            "{x:.8},{observation:.8},{:.8},{:.8},{:.8},{q05:.8},{q95:.8}",
            theta.mean, theta.sigma, theta.nu,
        );
    }
}

Интерфейсы различаются правилами заимствования. GAMLSS принимает &mut self, чтобы повторно использовать внутренние буферы, а функции обратного вызова argmin принимают &self. Локальный RefCell согласует эти интерфейсы без копирования рабочей области при каждом вызове.

В качестве алгоритма используется L-BFGS с памятью из десяти пар обновлений. MoreThuenteLineSearch подбирает длину шага с условиями Вольфе, поэтому шаг не задаётся вручную.

fn fit_once<O>(
    objective: O,
    constant_start: &[f64],
    nu_range: Range<usize>,
    nu_start: f64,
) -> Result<Fit, Error>
where
    O: Objective,
    O::Error: std::error::Error + Send + Sync + 'static,
{
    let mut initial_parameters = constant_start.to_vec();
    initial_parameters[nu_range].fill(nu_start);

    let problem = ArgminObjective::new(objective);
    let line_search = MoreThuenteLineSearch::new().with_c(1.0e-4, 0.9)?;
    let solver = LBFGS::new(line_search, LBFGS_MEMORY)
        .with_tolerance_grad(1.0e-5)?
        .with_tolerance_cost(1.0e-8)?;
    let result = Executor::new(problem, solver)
        .configure(|state| state.param(initial_parameters).max_iters(MAX_ITERATIONS))
        .run()?;

    let parameters = result
        .state
        .get_best_param()
        .or_else(|| result.state.get_param())
        .cloned()
        .ok_or_else(|| Error::msg("argmin finished without a parameter vector"))?;

    Ok(Fit {
        parameters,
        objective: result.state.get_best_cost(),
        iterations: result.state.get_iter(),
        termination: result
            .state
            .get_termination_reason()
            .map_or_else(|| "unknown".to_owned(), ToString::to_string),
        nu_start,
    })
}

Остановку задают допуски по норме градиента и изменению целевой функции, а также верхняя граница 750 итераций. После run() пример печатает причину остановки и заново вычисленную норму градиента.

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

cargo run --example a1_skew_normal

Для воспроизведения графика пример умеет печатать исходный отклик, оценённые параметры и границы 90%-интервала в формате CSV:

cargo run --quiet --example a1_skew_normal -- --csv > a1-fit.csv

Результат и диагностика

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

Интервал на рисунке вычислен по квантилям оценённого асимметричного нормального распределения. При $\nu (x)\ne 0$ он в общем случае не симметричен относительно $m(x)$, поэтому формула $m\pm 1,645s$ здесь неверна. Нижняя панель показывает, что направление асимметрии меняется по x: параметр $\nu$ положителен на краях диапазона и отрицателен в средней части.

    let fitted = model.predict_theta(&parameters)?;
    let intervals_90 = fitted
        .iter()
        .map(|theta| {
            (
                model.family().quantile(0.05, theta),
                model.family().quantile(0.95, theta),
            )
        })
        .collect::<Vec<_>>();
    if output_mode == OutputMode::Csv {
        print_prediction_csv(x, y, &fitted, &intervals_90);
        return Ok(());
    }

    let diagnostics = model.training_diagnostics(&parameters)?;
    let pit = model.pit_values(&parameters)?;
    let residuals = model.quantile_residuals(&parameters)?;
    let coverage_90 = intervals_90
        .iter()
        .zip(y)
        .filter(|(interval, observation)| (interval.0..=interval.1).contains(observation))
        .count() as f64
        / n as f64;

    let (mean_min, mean_max) = range(fitted.iter().map(|theta| theta.mean));
    let (sigma_min, sigma_max) = range(fitted.iter().map(|theta| theta.sigma));
    let (nu_min, nu_max) = range(fitted.iter().map(|theta| theta.nu));
    let (pit_mean, pit_sd) = mean_sd(&pit);
    let (residual_mean, residual_sd) = mean_sd(&residuals);

Среднее PIT для текущей модели близко к 0.5, а нормированные квантильные остатки имеют среднее около нуля и стандартное отклонение около единицы. Все приведённые значения вычислены на обучающих данных; в примере они служат компактной проверкой результата оптимизации.

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

Полный код примера находится в 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$. При большой размерности неструктурированная ковариационная модель быстро становится статистически и вычислительно дорогой; тогда нужны регуляризация, факторные или разреженные структуры и более сильные предположения о зависимости.

Двухкомпонентная гауссовская смесь для данных 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.

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

Литература и дополнительные материалы

Источники ниже полезны как статистический ориентир, коллекция удачных учебных ходов и набор пограничных случаев.

Книги и обзорные статьи

• Generalized Additive Models for Location, Scale and Shape: A Distributional Regression Approach, with Applications — современная книга 2024 года и открытый сайт-сопровождение с кодом, рисунками и наборами данных. Полезна разделением материала на основы, статистический вывод и полноценные прикладные разборы.
• Flexible Regression and Smoothing: Using GAMLSS in R — практическая книга о построении моделей, аддитивных членах, отборе, диагностике и оценивании центилей. Особенно полезен приём сравнивать несколько классов регрессии на одних данных.
• Distributions for Modeling Location, Scale, and Shape: Using GAMLSS in R — справочник по непрерывным, счётным и смешанным распределениям, их носителям, асимметрии и хвостам. Для этой книги он служит справочником, а не образцом линейного порядка глав.
• Rigby and Stasinopoulos (2005), Generalized additive models for location, scale and shape — исходная статья с общей постановкой GAMLSS, штрафным правдоподобием и аддитивными предикторами.
• Distributional regression using generalized additive models for location, scale and shape — недавний обзор 2026 года, связывающий GAMLSS с более широкой областью вероятностной регрессии.

Воспроизводимые примеры анализа

• gamlss2: First Steps моделирует распределение дневной максимальной температуры и превращает его в вероятность жаркого дня. Это хороший образец главы, которая заканчивается решением прикладного вопроса, а не таблицей коэффициентов.
• gamlss2: Centile (Quantile) Estimation проходит маршрут данные -> оценивание -> диагностика -> центильные кривые на dbbmi. Этот маршрут положен в основу запланированного прикладного разбора кривых роста.
• GAMLSS companion datasets публикует наборы данных из прикладных разборов в CSV для переноса в другие программы. Перед включением копии в репозиторий всё равно нужно проверить указание авторства и лицензию исходного набора.
• GAMLSS short course: Distributional regression показывает короткие сравнительные разборы и содержит полезные численные эталоны.

Диагностика и оценка прогнозов

• Dunn and Smyth (1996), Randomized quantile residuals объясняет квантильные остатки и необходимую рандомизацию для дискретных откликов. Это ключевой источник для будущего расширения диагностики счётных моделей.
• Gneiting and Raftery (2007), Strictly Proper Scoring Rules, Prediction, and Estimation даёт основу для логарифмической оценки, CRPS и корректного сравнения вероятностных прогнозов.

Наборы данных

• gamlss.data reference manual описывает dbbmi, abdom, speech и другие классические наборы данных вместе с происхождением и структурой переменных.
• UCI Bike Sharing содержит почасовые и посуточные числа поездок, погоду и календарные признаки и распространяется под CC BY 4.0.
• GasolineYield documentation описывает небольшой классический набор данных с долей выхода бензина строго внутри единичного интервала.
• French MTPL frequency analysis документирует задачу моделирования частоты страховых случаев для freMTPL/freMTPL2; таблицы частоты и тяжести убытка удобны для заключительного сквозного примера.