Надежность подобранной кривой?

Я хотел бы оценить неопределенность или надежность подобранной кривой. Я намеренно не называю точную математическую величину, которую ищу, потому что не знаю, что это.

Здесь (энергия) является зависимой переменной (отклик), а V (объем) является независимой переменной. Я хотел бы найти кривую энергии-объема, E ( V ) , некоторого материала. Поэтому я сделал несколько расчетов с компьютерной программой квантовой химии, чтобы получить энергию для некоторых объемов образцов (зеленые кружки на графике).$E$$V$$E(V)$

Затем я подобрал эти выборки данных с помощью функции Берча – Мурнагана : который зависит от четырех параметров: E 0 , V 0 , B 0 , B ′ 0 . Я также предполагаю, что это правильная функция подгонки, поэтому все ошибки происходят из-за шума выборок. В дальнейшем, насаженная функция ( Е ) будет записана как функция V .

$$\mathbb{E}(E|V) = E_0 + \frac{9V_0B_0}{16} \left\{ \left[\left(\frac{V_0}{V}\right)^\frac{2}{3}-1\right]^3B_0^\prime + \left[\left(\frac{V_0}{V}\right)^\frac{2}{3}-1\right]^2 \left[6-4\left(\frac{V_0}{V}\right)^\frac{2}{3}\right]\right\}\;,$$ $E_0, V_0, B_0, B_0'$$(\hat{E})$$V$

Здесь вы можете увидеть результат (в соответствии с алгоритмом наименьших квадратов). Переменной у-оси , а переменная х-оси V . Синяя линия - это подгонка, а зеленые кружки - точки отбора.$E$$V$

$\hat{E}(V)$

Мое интуиция говорит мне, что подобранная кривая наиболее надежна в середине, поэтому я предполагаю, что неопределенность (скажем, диапазон неопределенности) должна увеличиться ближе к концу выборочных данных, как на этом рисунке:

Однако, что это за мера, которую я ищу, и как я могу ее рассчитать?

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

Моя идея найти желаемую оценку неопределенности состоит в том, чтобы рассчитать следующую «ошибку» на основе параметров, которые вы изучаете в школе ( распространение неопределенности ):

$$\Delta E(V) = \sqrt{ \left(\frac{\partial E(V)}{\partial E_0} \Delta E_0\right)^2 + \left(\frac{\partial E(V)}{\partial V_0} \Delta V_0\right)^2 + \left(\frac{\partial E(V)}{\partial B_0} \Delta B_0\right)^2 + \left(\frac{\partial E(V)}{\partial B_0'} \Delta B_0'\right)^2}$$ $\Delta E_0, \Delta V_0, \Delta B_0$$\Delta B_0'$

Это приемлемый подход или я делаю это неправильно?

PS: я знаю, что я мог бы просто суммировать квадраты остатков между моими выборками данных и кривой, чтобы получить какую-то «стандартную ошибку», но это не зависит от объема.

Это обычная проблема наименьших квадратов!

определяющий

$$x = V^{-2/3}, \ w = V_0^{1/3},$$

модель можно переписать

$$\mathbb{E}(E|V) = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3$$

$\beta=(\beta_i)^\prime$

$$16 \beta = \pmatrix{16 E_0 + 54 B_0 w^3 - 9 B_0 B_0^\prime w^3\\ -144 B_0 w^5 + 27 B_0 B_0^\prime w^5\\ 126 B_0 w^7 - 27 B_0 B_0^\prime w^7\\ -36 B_0 w^9 + 9 B_0 B_0^\prime w^9}.$$

$B_0, B_0^\prime$$w$$B_0, B_0^\prime, w$$E_0$$\beta$

$(E_0, B_0, B_0^\prime, V_0)$$E$

$\hat\beta$R

#
# The data.
#
X <- data.frame(V=c(41, 43, 46, 48, 51, 53, 55.5, 58, 60, 62.5),
                E=c(-48.05, -48.5, -48.8, -49.03, -49.2, -49.3, -49.35, 
                    -49.34, -49.31, -49.27))
#
# OLS regression.
#
fit <- lm(E ~ I(V^(-2/3)) + I(V^(-4/3)) + I(V^(-6/3)), data=X)
summary(fit)
beta <- coef(fit)
#
# Prediction, including standard errors of prediction.
#
V0 <- seq(40, 65)
y <- predict(fit, se.fit=TRUE, newdata=data.frame(V=V0))
#
# Plot the data, the fit, and a three-SEP band.
#
plot(X$V, X$E, xlab="Volume", ylab="Energy", bty="n", xlim=c(40, 60))
polygon(c(V0, rev(V0)), c(y$fit + 3*y$se.fit, rev(y$fit - 3*y$se.fit)),
        border=NA, col="#f0f0f0")
curve(outer(x^(-2/3), 0:3, `^`) %*% beta, add=TRUE, col="Red", lwd=2)
points(X$V, X$E)

---

$\beta$

$I$$g$

$$-g^tIg$$ Это дает вам оценочную дисперсию для этой зависимой переменной. Возьмите квадратный корень, чтобы получить расчетное стандартное отклонение. тогда доверительные интервалы - это прогнозируемое значение + - два стандартных отклонения. Это стандартная вещь правдоподобия. для особого случая нелинейной регрессии вы можете скорректировать степени свободы. У вас есть 10 наблюдений и 4 параметра, так что вы можете увеличить оценку дисперсии в модели, умножив ее на 10/6. Несколько программных пакетов сделают это за вас. Я записал вашу модель в AD Model в AD Model Builder и подогнал ее и рассчитал (неизмененные) отклонения. Они будут немного отличаться от ваших, потому что мне пришлось немного угадать значения.

                    estimate   std dev
10   pred_E      -4.8495e+01 7.5100e-03
11   pred_E      -4.8810e+01 7.9983e-03
12   pred_E      -4.9028e+01 7.5675e-03
13   pred_E      -4.9224e+01 6.4801e-03
14   pred_E      -4.9303e+01 6.8034e-03
15   pred_E      -4.9328e+01 7.1726e-03
16   pred_E      -4.9329e+01 7.0249e-03
17   pred_E      -4.9297e+01 7.1977e-03
18   pred_E      -4.9252e+01 1.1615e-02

Это можно сделать для любой зависимой переменной в AD Model Builder. Каждый объявляет переменную в соответствующем месте в коде, как это

   sdreport_number dep

и пишет код для оценки зависимой переменной, как это

dep=sqrt(V0-cube(Bp0)/(1+2*max(V)));

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

19   dep          7.2535e+00 1.0980e-01

Я изменил программу, включив в нее код для расчета доверительных интервалов для функции объема энтальпии. Файл кода (TPL) выглядит следующим образом:

DATA_SECTION
 init_int nobs
 init_matrix data(1,nobs,1,2)
 vector E
 vector V
 number Vmean
LOC_CALCS
 E=column(data,2);
 V=column(data,1);
 Vmean=mean(V);

PARAMETER_SECTION
 init_number E0
 init_number log_V0_coff(2)
 init_number log_B0(3)
 init_number log_Bp0(3)
 init_bounded_number a(.9,1.1)
 sdreport_number V0
 sdreport_number B0
 sdreport_number Bp0
 sdreport_vector pred_E(1,nobs)
 sdreport_vector P(1,nobs)
 sdreport_vector H(1,nobs)
 sdreport_number dep
 objective_function_value f
PROCEDURE_SECTION
  V0=exp(log_V0_coff)*Vmean;
  B0=exp(log_B0);
  Bp0=exp(log_Bp0);
  if (current_phase()<4)
  f+=square(log_V0_coff) +square(log_B0);

  dvar_vector sv=pow(V0/V,0.66666667);
  pred_E=E0 + 9*V0*B0*(cube(sv-1.0)*Bp0
    + elem_prod(square(sv-1.0),(6-4*sv)));

  dvar_vector r2=square(E-pred_E);
  dvariable vhat=sum(r2)/nobs;
  dvariable v=a*vhat;
  f=0.5*nobs*log(v)+sum(r2)/(2.0*v);

  // code to calculate the  enthalpy-volume function
  double delta=1.e-4;
  dvar_vector svp=pow(V0/(V+delta),0.66666667);
  dvar_vector svm=pow(V0/(V-delta),0.66666667);
  P = -((9*V0*B0*(cube(svp-1.0)*Bp0
      + elem_prod(square(svp-1.0),(6-4*svp))))
      -(9*V0*B0*(cube(svm-1.0)*Bp0
      + elem_prod(square(svm-1.0),(6-4*svm)))))/(2.0*delta);
  H=E+elem_prod(P,V);

dep=sqrt(V0-cube(Bp0)/(1+2*max(V)));

Затем я переоборудовал модель, чтобы получить стандартные разработки для оценок H.

29   H           -3.9550e+01 5.9163e-01
30   H           -4.1554e+01 2.8707e-01
31   H           -4.3844e+01 1.2333e-01
32   H           -4.5212e+01 1.5011e-01
33   H           -4.6859e+01 1.5434e-01
34   H           -4.7813e+01 1.2679e-01
35   H           -4.8808e+01 1.1036e-01
36   H           -4.9626e+01 1.8374e-01
37   H           -5.0186e+01 2.8421e-01
38   H           -5.0806e+01 4.3179e-01

Они рассчитаны для ваших наблюдаемых значений V, но могут быть легко рассчитаны для любого значения V.

Было отмечено, что на самом деле это линейная модель, для которой существует простой R-код для оценки параметров через OLS. Это очень привлекательно, особенно для наивных пользователей. Однако, поскольку работа Хубера более тридцати лет назад, мы знаем или должны знать, что, вероятно, почти всегда следует заменить OLS на умеренно надежную альтернативу. Я полагаю, что причина этого обычно не в том, что надежные методы по своей природе нелинейны. С этой точки зрения простые привлекательные методы OLS в R являются скорее ловушкой, а не особенностью. Преимуществом подхода AD Model Builder является встроенная поддержка нелинейного моделирования. Чтобы изменить код наименьших квадратов на устойчивую нормальную смесь, нужно изменить только одну строку кода. Линия

    f=0.5*nobs*log(v)+sum(r2)/(2.0*v);

изменено на

f=0.5*nobs*log(v)
  -sum(log(0.95*exp(-0.5*r2/v) + 0.05/3.0*exp(-0.5*r2/(9.0*v))));

Степень избыточной дисперсии в моделях измеряется параметром a. Если a равно 1,0, дисперсия такая же, как для нормальной модели. Если имеется отклонение от инфляции по выбросам, мы ожидаем, что a будет меньше 1,0. Для этих данных оценка a составляет около 0,23, так что дисперсия составляет около 1/4 дисперсии для нормальной модели. Интерпретация заключается в том, что выбросы увеличили оценку дисперсии примерно в 4 раза. Результатом этого является увеличение размера доверительных границ для параметров для модели OLS. Это представляет собой потерю эффективности. Для модели нормальной смеси оцененные стандартные отклонения для функции энтальпии-объема

 29   H           -3.9777e+01 3.3845e-01
 30   H           -4.1566e+01 1.6179e-01
 31   H           -4.3688e+01 7.6799e-02
 32   H           -4.5018e+01 9.4855e-02
 33   H           -4.6684e+01 9.5829e-02
 34   H           -4.7688e+01 7.7409e-02
 35   H           -4.8772e+01 6.2781e-02
 36   H           -4.9702e+01 1.0411e-01
 37   H           -5.0362e+01 1.6380e-01
 38   H           -5.1114e+01 2.5164e-01

Видно, что в точечных оценках есть небольшие изменения, в то время как доверительные пределы были уменьшены примерно до 60% от тех, которые были получены с помощью OLS.

Главное, что я хочу подчеркнуть, это то, что все измененные вычисления выполняются автоматически, как только одна строка кода изменяется в файле TPL.

Перекрестная проверка - это простой способ оценить надежность вашей кривой: https://en.wikipedia.org/wiki/Cross-validation_(statistics)

$\Delta E_{0},\Delta V_{0},\Delta B_{0}$$\Delta B'$

Вы можете рассчитать 1-кратную ошибку проверки, оставив одну из ваших точек вне аппроксимации и используя подобранную кривую для прогнозирования значения точки, которая была оставлена. Повторите это для всех точек так, чтобы каждый остался один раз. Затем вычислите ошибку проверки вашей окончательной кривой (кривой, снабженной всеми точками) как среднее значение ошибок прогнозирования.

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

Кроме того, вы можете отобразить ошибки прогнозирования в зависимости от объема, если хотите.