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

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

Сплайн-модель асимметричного нормального распределения для a1

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

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

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

$Y_{i} ∣ x_{i} \sim SN_{mean/sd} (m_{i}, s_{i}, ν_{i}), E [Y_{i} ∣ x_{i}] = m_{i}, SD (Y_{i} ∣ x_{i}) = s_{i} .$

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

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

$m_{i} = η_{m, i}, s_{i} = exp (clamp (η_{s, i}, - 8, 2)), ν_{i} = η_{ν, i} .$

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

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

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

$η_{r} (x) = j = 0 \sum K_{r} - 1 β_{r, j} B_{r, j} (x), r \in {m, s, ν} .$

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

Параметр	Маркер	Число базисных функций	Link-функция	Вес штрафа
Условное среднее $m (x)$	`Mean`	20	`Identity`	0.05
Стандартное отклонение $s (x)$	`Sigma`	14	`ClampedLog<-8, 2>`	0.10
Асимметрия $ν (x)$	`Nu`	12	`Identity`	0.01

Всего оптимизатор видит $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} (β) = \frac{λ}{K - 2} j = 0 \sum K - 3 (β_{j + 2} - 2 β_{j + 1} + β_{j})^{2} .$

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

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

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

$L (β) = \frac{1}{n} i = 1 \sum n - lo g p (y_{i} ∣ m_{i}, s_{i}, ν_{i}) + J_{m} + J_{s} + J_{ν} .$

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

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

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

Ненулевой старт здесь принципиален. В параметризации через среднее и стандартное отклонение $ν = 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

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

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

Величина	Результат
Усреднённая целевая функция со штрафом	-0.7715
Среднее отрицательное логарифмическое правдоподобие без штрафов	-0.7931
Сумма штрафов	0.0215
Норма градиента	$1.6 \cdot 1 0^{- 4}$
Диапазон $m (x)$	$[- 0.478, 0.917]$
Диапазон $s (x)$	$[0.050, 0.175]$
Диапазон $ν (x)$	$[- 2.91, 7.33]$
Покрытие 90%-интервала на обучающих данных	0.888

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

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

Keyboard shortcuts

GAMLSS на Rust