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

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