Степенная регрессия и Box–Cox: почему важны p0 и якобиан

Разбираем нелинейную регрессию: степенная модель и форма Box–Cox. Почему без p0 и якобиана оптимизация нестабильна и как влияет масштабирование на результат.

Нелинейная аппроксимация кривых может быть коварной, когда алгебраически эквивалентные модели ведут себя по‑разному под управлением оптимизатора. Типичный пример — регрессия степенного вида y = a + b * x ** c и её перепараметризация в духе Box–Cox. Вопрос в том, какая запись даст лучший результат. Короткий ответ с практической точки зрения оптимизации: без хорошего начального приближения и якобиана вы проверяете удачу оптимизатора, а не качество модели. Если задать и то и другое, эквивалентные формулировки сходятся одинаково.

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

Сравниваются две математически эквивалентные параметризации. Первая: y = a + b * x ** c. Вторая: d + e * (x ** f - 1) / f, при взаимосвязи параметров a = d - e / f, b = e / f, c = f. Было замечено, что первый вариант превышал значение maxfev по умолчанию и чувствителен к масштабированию признака, тогда как второй сходился с настройками по умолчанию, давал лучшее соответствие и почти не реагировал на масштаб. После масштабирования x до единичного стандартного отклонения обе записи дали одинаковое решение.

Воспроизводимый пример, показывающий расхождение

Следующий код иллюстрирует это поведение: обе формы подгоняются дважды — на исходном x и на масштабированном — с последующим сравнением MAE. Логика полностью повторяет описанную постановку.

from io import StringIO
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.optimize import curve_fit
from sklearn.metrics import mean_absolute_error
from sklearn.preprocessing import StandardScaler
def base_power(x, p, q, r):
    return p + q * x ** r
def shifted_boxcox(x, u, v, w):
    return u + v * (x ** w - 1) / w
DATA_SRC = StringIO(
    """x,y
410.0,73.06578085929756
25.0,205.29417389522575
72.0,110.48653325137172
51.0,168.52111516008628
15.0,119.75684720989004
72.0,164.46280991735537
73.0,145.53126751391522
36.0,161.41429319490925
40.0,219.91735537190084
190.0,89.18717804003897
21.0,203.3350571291969
47.0,170.12964176670877
21.0,197.1020714822368
18.0,96.43526170798899
53.0,117.55060034305319
50.0,189.89358650365415
43.0,179.03132807995385
863.0,69.63888057656149
71.0,131.42730764753813
40.0,205.2892561983471
65.0,131.3857292219426
401.0,133.81511047189076
50.0,115.65603442387814
58.0,151.99074870050802
50.0,165.8640803223824
21.0,210.87942861045792
236.0,124.21734182739671
53.0,180.11451429366744
12.0,320.77043917765184
36.0,244.3526170798898
25.0,202.41568198893515
21.0,184.03895128162597
29.0,165.64724945771087
25.0,218.1818181818182
72.0,161.8457300275482
130.0,107.38232466256511
84.0,177.52397865095088
38.0,57.524112172378224
50.0,168.132777815723
25.0,202.41568198893515
21.0,244.3978260657449
48.0,168.3167133392528
200.0,122.8403797554812
37.0,167.84185838731295
83.0,173.75445583988846
13.0,315.835929660122
11.0,314.47181327653976
32.0,203.68741889215215
200.0,123.96694214876034
39.0,110.59353869271226
39.0,190.81504521023686
40.0,235.53719008264466
37.0,181.71111484409758
25.0,215.55576804863057
40.0,235.53719008264466"""
)
frame = pd.read_csv(DATA_SRC)
scaler = StandardScaler(with_mean=False)
frame["x_std"] = scaler.fit_transform(frame[["x"]])
fig, axes = plt.subplots(1, 2, figsize=(10, 6), sharey=True)
for idx, use_scaled in enumerate([False, True]):
    feat_col = "x"
    if use_scaled:
        feat_col = "x_std"
    axes[idx].scatter(frame[feat_col], frame["y"], label="Data", alpha=0.5)
    for variant_idx, (fn, style) in enumerate(
        zip([base_power, shifted_boxcox], ["--", ":"])
    ):
        theta, cov = curve_fit(
            fn,
            frame[feat_col],
            frame["y"],
            maxfev=10_000,
        )
        frame["pred"] = fn(frame[feat_col], *theta)
        mae_val = mean_absolute_error(frame["y"], frame["pred"])
        grid = np.linspace(frame[feat_col].min(), frame[feat_col].max(), 100)
        yhat_line = fn(grid, *theta)
        axes[idx].plot(
            grid,
            yhat_line,
            linestyle=style,
            label=f"Form {variant_idx + 1} (MAE: {mae_val:.2f})",
        )
    axes[idx].legend()
    axes[idx].set_title(f"Scaled: {use_scaled}")
    axes[idx].set_xlabel(feat_col)
    axes[idx].set_ylabel("y")
plt.show()

Почему эквивалентные модели здесь расходятся

Корневая причина — отсутствие начальной оценки p0. Без неё оптимизатор стартует с внутренних значений по умолчанию. Эквивалентные функции способны уводить решатель в разные области параметрического пространства, обеспечивать различную скорость сходимости и приводить к разным локальным минимумам. Это объясняет, почему одна форма будто бы требовала большего maxfev и реагировала на масштабирование, а другая — нет. Когда x стандартизирован до единичной дисперсии, оба варианта совпадают, потому что такое масштабирование фактически выравнивает области поиска по параметрам.

Как заставить оптимизатор работать предсказуемо: задайте начальную оценку

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

import io
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
def base_curve(x: np.ndarray, aa: float, bb: float, cc: float) -> np.ndarray:
    return aa + bb * x ** cc
def run_fit(x: np.ndarray, y: np.ndarray, guess: np.typing.ArrayLike) -> np.ndarray:
    popt, pcov, info, msg, status = curve_fit(
        f=base_curve, xdata=x, ydata=y, p0=guess, full_output=True,
    )
    if status not in {1, 2, 3, 4}:
        raise ValueError(msg)
    print(msg, f'in {info["nfev"]} calls')
    return popt
def draw_fig(
    x: np.ndarray, y: np.ndarray,
    theta_guess: np.ndarray,
    theta_fit: np.ndarray,
) -> plt.Figure:
    xgrid = np.logspace(start=1, stop=3, num=200)
    fig, ax = plt.subplots()
    ax.set_xscale('log')
    ax.scatter(x, y, label='Data')
    ax.plot(xgrid, base_curve(xgrid, *theta_guess), label='Est')
    ax.plot(xgrid, base_curve(xgrid, *theta_fit), label='Fit')
    ax.legend()
    return fig
def run_demo() -> None:
    raw = io.StringIO(
'''410.0,73.06578085929756
25.0,205.29417389522575
72.0,110.48653325137172
51.0,168.52111516008628
15.0,119.75684720989004
72.0,164.46280991735537
73.0,145.53126751391522
36.0,161.41429319490925
40.0,219.91735537190084
190.0,89.18717804003897
21.0,203.3350571291969
47.0,170.12964176670877
21.0,197.1020714822368
18.0,96.43526170798899
53.0,117.55060034305319
50.0,189.89358650365415
43.0,179.03132807995385
863.0,69.63888057656149
71.0,131.42730764753813
40.0,205.2892561983471
65.0,131.3857292219426
401.0,133.81511047189076
50.0,115.65603442387814
58.0,151.99074870050802
50.0,165.8640803223824
21.0,210.87942861045792
236.0,124.21734182739671
53.0,180.11451429366744
12.0,320.77043917765184
36.0,244.3526170798898
25.0,202.41568198893515
21.0,184.03895128162597
29.0,165.64724945771087
25.0,218.1818181818182
72.0,161.8457300275482
130.0,107.38232466256511
84.0,177.52397865095088
38.0,57.524112172378224
50.0,168.132777815723
25.0,202.41568198893515
21.0,244.3978260657449
48.0,168.3167133392528
200.0,122.8403797554812
37.0,167.84185838731295
83.0,173.75445583988846
13.0,315.835929660122
11.0,314.47181327653976
32.0,203.68741889215215
200.0,123.96694214876034
39.0,110.59353869271226
39.0,190.81504521023686
40.0,235.53719008264466
37.0,181.71111484409758
25.0,215.55576804863057
40.0,235.53719008264466''')
    with raw:
        x, y = np.loadtxt(raw, delimiter=',').T
    theta0 = (0, 1000, -0.5)
    theta = run_fit(x, y, theta0)
    print('theta =', theta)
    draw_fig(x, y, theta_guess=theta0, theta_fit=theta)
    plt.show()
if __name__ == '__main__':
    run_demo()

И фактическое, и прогнозируемое относительное уменьшение суммы квадратов не превышают 0.000000 за 25 вызовов

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

Шаг дальше: передайте якобиан

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

import io
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
def base_curve(x: np.ndarray, aa: float, bb: float, cc: float) -> np.ndarray:
    return aa + bb * x ** cc
def grad_fn(x: np.ndarray, aa: float, bb: float, cc: float) -> np.ndarray:
    x_pow = x ** cc
    j = np.empty(shape=(3, x.size), dtype=x.dtype)
    j[0] = 1
    j[1] = x_pow
    j[2] = bb * x_pow * np.log(x)
    return j
def run_fit_with_jac(x: np.ndarray, y: np.ndarray, guess: np.typing.ArrayLike) -> np.ndarray:
    popt, pcov, info, msg, status = curve_fit(
        f=base_curve, jac=grad_fn, xdata=x, ydata=y, p0=guess,
        full_output=True, col_deriv=True,
    )
    print(f'status={status}, {info["nfev"]} model calls, {info["njev"]} Jacobian calls')
    print(msg)
    return popt
def draw_fig(
    x: np.ndarray, y: np.ndarray,
    theta_guess: np.typing.ArrayLike,
    theta_fit: np.typing.ArrayLike,
) -> plt.Figure:
    xgrid = np.logspace(start=1, stop=3, num=100)
    fig, ax = plt.subplots()
    ax.set_xscale('log')
    ax.scatter(x, y, alpha=0.5, label='Data')
    ax.plot(xgrid, base_curve(xgrid, *theta_guess), label='Est')
    ax.plot(xgrid, base_curve(xgrid, *theta_fit), label='Fit')
    ax.legend()
    return fig
def run_demo_jac() -> None:
    raw = io.StringIO(
'''410.0,73.06578085929756
25.0,205.29417389522575
72.0,110.48653325137172
51.0,168.52111516008628
15.0,119.75684720989004
72.0,164.46280991735537
73.0,145.53126751391522
36.0,161.41429319490925
40.0,219.91735537190084
190.0,89.18717804003897
21.0,203.3350571291969
47.0,170.12964176670877
21.0,197.1020714822368
18.0,96.43526170798899
53.0,117.55060034305319
50.0,189.89358650365415
43.0,179.03132807995385
863.0,69.63888057656149
71.0,131.42730764753813
40.0,205.2892561983471
65.0,131.3857292219426
401.0,133.81511047189076
50.0,115.65603442387814
58.0,151.99074870050802
50.0,165.8640803223824
21.0,210.87942861045792
236.0,124.21734182739671
53.0,180.11451429366744
12.0,320.77043917765184
36.0,244.3526170798898
25.0,202.41568198893515
21.0,184.03895128162597
29.0,165.64724945771087
25.0,218.1818181818182
72.0,161.8457300275482
130.0,107.38232466256511
84.0,177.52397865095088
38.0,57.524112172378224
50.0,168.132777815723
25.0,202.41568198893515
21.0,244.3978260657449
48.0,168.3167133392528
200.0,122.8403797554812
37.0,167.84185838731295
83.0,173.75445583988846
13.0,315.835929660122
11.0,314.47181327653976
32.0,203.68741889215215
200.0,123.96694214876034
39.0,110.59353869271226
39.0,190.81504521023686
40.0,235.53719008264466
37.0,181.71111484409758
25.0,215.55576804863057
40.0,235.53719008264466''')
    with raw:
        x, y = np.loadtxt(raw, delimiter=',').T
    theta0 = 0, 1000, -0.5
    aa, bb, cc = theta = run_fit_with_jac(x, y, theta0)
    print(f'MAE = {np.abs(y - base_curve(x, *theta)).mean():.2f}')
    print(f'a = {aa:+.3f}')
    print(f'b = {bb:+.2f}')
    print(f'c = {cc:+.4f}')
    print(f'd = {aa+bb:+.2f}')
    print(f'e = {bb*cc:+.2f}')
    print(f'f = {cc:+.4f}')
    draw_fig(x, y, theta_guess=theta0, theta_fit=theta)
    plt.show()
if __name__ == '__main__':
    run_demo_jac()

status=1, 7 вызовов модели, 6 вызовов якобиана. И фактическое, и прогнозируемое относительное уменьшение суммы квадратов не превышают 0.000000. MAE = 29.09.

Это показывает, что при наличии начальной оценки и якобиана подгонка для базовой степенной модели проходит быстро и устойчиво. По отображению параметров a = d − e / f, b = e / f, c = f форма в стиле Box–Cox при эквивалентных начальных значениях выйдет на то же решение.

Почему это важно

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

Ещё одно наблюдение: добавление одной степени свободы в виде a + b * (x − c) ** d снизило MAE с 29.1 до 27.9 на тех же данных. Насколько такая модель уместна — зависит от задачи, и сравнение с двумя эквивалентными параметризациями выше нельзя считать полностью корректным.

Выводы и рекомендации

При подгонке y = a + b * x ** c или её перепараметризации в стиле Box–Cox в первую очередь задавайте продуманные начальные оценки и, по возможности, передавайте якобиан. При этих условиях алгебраически эквивалентные записи не будут заметно отличаться по качеству, сходимости и чувствительности к масштабированию. Более гибкие модели стоит рассматривать лишь тогда, когда они отражают реальный процесс, и сравнивать их на одинаковых условиях.