Подгонка болцмановской модели в SciPy: скобки, масштаб и kB

Почему в SciPy нелинейная подгонка болцмановской модели даёт «covariance could not be estimated», и как это починить: скобки, масштаб kB, свободная амплитуда N.

Подгонка нелинейных моделей к экспериментальным кривым — привычная задача, но пара тонких ловушек способна заставить оптимизатор застрять или вернуть бессмысленные параметры. Ниже — компактный разбор реального случая: одномерная подгонка по выражению болцмановского типа приводит к сообщению «covariance could not be estimated», сходится к T=1 и даёт вводящий в заблуждение график. Разберёмся, что пошло не так, и как это поправить, не меняя задуманной физики.

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

Целевая модель описывается выражением, пропорциональным B(2J+1)/(kB·T)·exp(−B·J(J+1)/(kB·T)). Независимая переменная — J, неизвестный параметр — T, а kB и B — константы. Прямой заход на подгонку в SciPy может выглядеть так:

# Константы
k_b_const = 1.380649e-23
beta_coeff = 0.3808
# Функция подгонки
def fit_fun(j_var, t_param):
    return (beta_coeff * (2 * j_var + 1) / k_b_const * t_param *
            np.exp(-beta_coeff * j_var * (j_var + 1) / k_b_const * t_param))
p_opt, p_cov = optimize.curve_fit(fit_fun, x_vals, y_vals)
print(p_opt)
plt.plot(x_vals, y_vals, "o")
plt.plot(x_vals, fit_fun(x_vals, p_opt[0]))

Оптимизатор отвечает ошибкой про ковариацию параметров, а единственный параметр становится равным 1:

OptimizeWarning: Covariance of the parameters could not be estimated

Вот набор данных, использованный для подгонки:

x_vals = [2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48]
y_vals = [0.185,0.375,0.548,0.695,0.849,0.931,0.996,0.992,1.0,0.977,0.927,0.834,0.691,0.68,0.575,0.479,0.421,0.351,0.259,0.208,0.162,0.135,0.093,0.066]

Что на самом деле идёт не так

Первая проблема — приоритет операций. Запись a/kb*T — это не то же самое, что a/(kb*T). В исходной форме и множитель, и экспонента фактически умножаются на T, а не делятся на (kB·T). В такой конфигурации оптимизатор по сути учит 1/T вместо T. Это само по себе не объясняет сбой ковариации, но искажает параметризацию и усиливает нестабильность.

Вторая проблема — нормировка целевого вектора. Переданные y-значения образуют распределение с единичной площадью под кривой. Если в модели жёстко зафиксировать амплитуду равной единице, кривая не сможет подстроиться под наблюдаемый масштаб, и сходимость блокируется. Добавление явного масштабирующего параметра N решает это: модель подбирает общий уровень, а T задаёт форму распределения.

Третья проблема — согласованность единиц и масштабирование. Константа Больцмана задана в Дж/К, тогда как единицы измерения B не указаны. Показатель экспоненты должен быть безразмерным. Использование kB в эВ/К вместо Дж/К даёт более подходящий численный масштаб в данном контексте и обеспечивает стабильные, правдоподобные оценки T.

Есть и две практические мелочи. Всегда ставьте скобки вокруг знаменателей вроде (kB*T), чтобы избежать ошибок приоритета. И при построении графика подставляйте тот же массив x, что использовался при подгонке: обращение к неопределённой или иной переменной j ведёт к несоответствующим или вводящим в заблуждение графикам.

Пошаговое исправление и рабочий пример

Лекарство минимально: добавить скобки вокруг kB*T, ввести свободную амплитуду N и перейти на kB в эВ/К. Всё остальное остаётся прежним. Ниже — полный самодостаточный пример, который подгоняет предоставленные данные и валидирует результат.

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize, integrate
from sklearn import metrics
# kB в эВ/К для корректного масштабирования
kB_ev = 8.617333262e-5
B_const = 0.3808
def shape_fn(j_idx, temp, amp):
    return amp * (
        B_const * (2 * j_idx + 1) / (kB_ev * temp) *
        np.exp(- B_const * j_idx * (j_idx + 1) / (kB_ev * temp))
    )
j_data = np.array([2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48])
y_data = np.array([0.185,0.375,0.548,0.695,0.849,0.931,0.996,0.992,1.0,0.977,0.927,0.834,0.691,0.68,0.575,0.479,0.421,0.351,0.259,0.208,0.162,0.135,0.093,0.066])
# Начальное приближение [T, N]
params_hat, cov_hat = optimize.curve_fit(shape_fn, j_data, y_data, p0=[1e5, 1])
# Неопределенности из ковариационной матрицы
stderr = np.sqrt(np.diag(cov_hat))
# Прогноз и R^2
pred = shape_fn(j_data, *params_hat)
r2 = metrics.r2_score(y_data, pred)
# Гладкая кривая для визуализации
j_grid = np.linspace(j_data.min(), j_data.max(), 200)
y_grid = shape_fn(j_grid, *params_hat)
fig, ax = plt.subplots()
ax.scatter(j_data, y_data)
ax.plot(j_grid, y_grid)
ax.set_xlabel("J'')")
ax.set_ylabel("N(J'')")
ax.grid()
print("Estimated [T, N]:", params_hat)
print("Std errors        :", stderr)
print("R2                :", r2)

В этой конфигурации регрессия сходится корректно. Получаются примерно T = 2.51e6 ± 0.03e6 и N = 2.76e1 ± 0.02e1, а качество аппроксимации высокое: R² ≈ 0.994.

Если же использовать kB в Дж/К, оценка температуры окажется порядка 1.6e+25 — в данном месте это численно неудобно. Причина исключительно в единицах и масштабировании: показатель экспоненты должен оставаться безразмерным.

Поскольку данные представляют нормированное распределение, можно проверить, что подогнанная кривая интегрируется к единице после удаления амплитуды. Это легко проверить численно:

j_check = np.linspace(0, 100, 2000)
y_check = shape_fn(j_check, *params_hat)
area = integrate.trapezoid(y_check / params_hat[-1], j_check)
print("Area under curve (normalized):", area)

Площадь фактически равна единице, что согласуется с интерпретацией y как распределения с единичной площадью.

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

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

Практические выводы

Ставьте скобки в знаменателях вроде (kB*T) и в множителе, и в экспоненте. Добавляйте свободную амплитуду N при подгонке распределений, у которых площадь фиксирована конструктивно. Держите показатель экспоненты безразмерным, согласуя единицы; использование kB в эВ/К в данном случае — прагматичный выбор, дающий устойчивые оценки. И, наконец, будьте последовательны при визуализации: вычисляйте подогнанную функцию на том же массиве x, что применялся при подгонке, чтобы избежать несоответствий.

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