На одном из предыдущих занятий было рассмотрено получение параметров нединейной функции, фитующей экпериментальнеы данные в смысле наименьших квадратов. Однако, использованное там примение функции k = lsqcurvefit(demodel,k0,SOC1dn(ind),Vd1(ind)); возвращало только сами значения параметров k, в то время как более статистически-убедительное их представление должно сопровождаться доверительным интервалов при заданном уровне достоверности. Помимо этого, имеет смысл также указывать и интервал, в который могут попасть предсказания использованной модели, принимая во внимание данные доверительные интервалы параметров.
В GNU Octave (это же работает и в MATLAB, хотя там можно использовать и готовые функции fit и predint) для этого следует применять расширенный надор выходных параметров функции lsqcurvefit (здесь и далее будут приведены долько дополнения к предыдущему коду, все остальные имена переменных и промежудточные вычисления - те же.):
[k, resnorm, residual, exitflag, output, lambda,jacobian] = lsqcurvefit(demodel,k0,SOC1dn(ind),Vd1(ind),[],[]);
Здесь существенно, что в аргумене функции необходими как минимум добавить еще два входных аргумента - верхние и нижние границы параметров; если нет необходимости их задвать конкретными числами, то можно просто поставить пустые переменные [], как сделано выше.
Расширенный набор выходных параметров (их матемаматичсекий смысл см. по ссылке) используется для расчета следующих вспомогательных величин:
Количество точек фитуемых данных:
n = length(ind);
Число параметров модели:
p=length(k);
Следующее из этой пары количество степеней свободы (что существенно - число параметров должно быть больше, чем число данных, иначе при их совпадении модельная кривая просто проходит по всем точкам, что убивает всякую неопредленность и статистически малоосмысленно, а пареметров больше, то модель не может функционировать как таковая):
df = n-p;
Статистическая значимость, выраженная числом от 0 до 1 (т.е., например, 0.05 для 95% доверительных интервалов)
alpha = 0.05;
Критическое значение распределения Стьюдента при данных статистической значимости и числе степеней свободы:
t_crit = tinv(1 - alpha/2, df);
Среднеквадратичное отклонение:
MSE = resnorm / (n - p);
Ковариационная матрица параметров:
CovB = MSE * inv(jacobian' * jacobian);
SE = sqrt(diag(CovB))';
Значения этих величин дают искомые нижную и верхнюю границы доверительного интервала для найденных параметров при заданном уровне надежности:
CI_lower = k - t_crit * SE;
CI_upper = k + t_crit * SE;
В частности, если сами параметры, найденные в смысле наилучшего приближения в среднеквадратичном смысле равны
,
до
Кроме доверительных интервалов для параметров, можно рассчитать и доверительные интервалы для точек, рассчитываемых по использованной модели (здесь расчет для тех же точек, в которых изерялись экпериментальные данные):
delta_conf = zeros(size(ind));
for i = 1:length(ind)
% Ji - якобиан для i-й точки наблюдения
Ji = jacobian(i, :);
se_fit = sqrt(MSE * (Ji * CovB * Ji'));
delta_conf(i) = t_crit * se_fit;
end
Более надежная велична - предсказательный интервал, вычисляется образом, похожим на расчет интервала параметров, т.е. с учетом дополнительного разброса, учитываемого распределением Стьюдента:
delta_pred = delta_conf + t_crit * sqrt(MSE);
При построении графика он задает коридор линий выше и ниже модельной кривой:
plot(SOC1dn(ind),demodel(k,SOC1dn(ind)),'-','LineWidth',1)
hold on
plot(SOC1dn(ind),demodel(k,SOC1dn(ind))- delta_pred,'--','LineWidth',1,'color','red')
plot(SOC1dn(ind),demodel(k,SOC1dn(ind))+ delta_pred,'--','LineWidth',1,'color','red')
Результрующий график:
На нижней панели представлен увеличенный кусок графика, который показывает, что экспериментальные точки действительно попадаеют в предстказательный интервал данной модели на 95% уровне статистической значимости.
Более подробное обсуждение сущности доверительных интервалов при нелинейной регресии, а также другие возможные варианты их определения (более подходящие, например, в случае, когда интервал несимметричен), можно найти в статье:
- O’Brien T. E., Silcox J. W. Nonlinear Regression Modelling: A Primer with Applications and Caveats. Bulletin of Mathematical Biology. 2024. 86 (2024) 40 [ссылка].
