¿Fiabilidad de una curva ajustada?

Question

Me gustaría estimar la incertidumbre o la fiabilidad de una curva ajustada. Intencionalmente no nombro una cantidad matemática precisa que estoy buscando, ya que no sé lo que es.

Aquí $E$ (energía) es la variable dependiente (respuesta) y $V$ (volumen) es la variable independiente. Me gustaría encontrar la curva Energía-Volumen, $E(V)$ de algún material. Así que hice algunos cálculos con un programa informático de química cuántica para obtener la energía de algunos volúmenes de muestra (círculos verdes en el gráfico).

A continuación, ajusté estas muestras de datos con el Función Birch-Murnaghan : $$ \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\}\;, $$ que depende de cuatro parámetros: $E_0, V_0, B_0, B_0'$ . También asumo que esta es la función de ajuste correcta, por lo que todos los errores sólo provienen del ruido de las muestras. En lo que sigue, la función ajustada $(\hat{E})$ se escribirá en función de $V$ .

Aquí puedes ver el resultado (ajuste con un algoritmo de mínimos cuadrados). La variable del eje Y es $E$ y la variable del eje x es $V$ . La línea azul es el ajuste y los círculos verdes son los puntos de la muestra.

Ahora necesito alguna medida de la fiabilidad (en el mejor de los casos en función del volumen) de esta curva ajustada, $\hat{E}(V)$ porque lo necesito para calcular otras magnitudes como las presiones de transición o las entalpías.

Mi intuición me dice que la curva ajustada es más fiable en el centro, así que supongo que la incertidumbre (digamos el rango de incertidumbre) debería aumentar cerca del final de los datos de la muestra, como en este esquema:

Sin embargo, ¿qué es este tipo de medida que busco y cómo puedo calcularla?

Para ser precisos, aquí sólo hay una fuente de error: Las muestras calculadas tienen ruido debido a los límites computacionales. Así que si calculara un conjunto denso de muestras de datos, formarían una curva irregular.

Mi idea para encontrar la estimación de la incertidumbre deseada es calcular el siguiente ''error'' basado en los parámetros como se aprende en la escuela ( propagación de la incertidumbre ):

$$ \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} $$ El $\Delta E_0, \Delta V_0, \Delta B_0$ y $\Delta B_0'$ son dadas por el software de adaptación.

¿Es un enfoque aceptable o lo estoy haciendo mal?

PD: Sé que también podría sumar los cuadrados de los residuos entre mis muestras de datos y la curva para obtener algún tipo de ''error estándar'', pero esto no depende del volumen.

Answer 1

Se trata de un problema ordinario de mínimos cuadrados.

Definición de

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

el modelo se puede reescribir

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

donde los coeficientes $\beta=(\beta_i)^\prime$ están relacionados algebraicamente con los coeficientes originales a través de

$$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}.$$

Esto es sencillo de resolver algebraicamente o numéricamente: elige la solución para la que $B_0, B_0^\prime$ y $w$ son positivos. La única razón para hacerlo es obtener estimaciones de $B_0, B_0^\prime, w$ y $E_0$ y para verificar que son físicamente significativos. Todos los análisis del ajuste pueden llevarse a cabo en términos de $\beta$ .

Este enfoque no sólo es mucho más sencillo que el ajuste no lineal, sino que también es más preciso: la matriz de varianza-covarianza para $(E_0, B_0, B_0^\prime, V_0)$ devuelto por un ajuste no lineal es sólo una aproximación cuadrática local a la distribución de muestreo de estos parámetros, mientras que (para los errores normalmente distribuidos en la medición $E$ En cualquier caso, los resultados de OLS no son aproximaciones.

Los intervalos de confianza, los intervalos de predicción, etc. pueden obtenerse de la forma habitual sin necesidad de encontrar estos valores: calcularlos en función de las estimaciones $\hat\beta$ y su matriz de varianza-covarianza. (¡Incluso Excel podría hacer esto!) Aquí hay un ejemplo, seguido de la (simple) R código que lo produjo.

#
# 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)

---

Si está interesado en la distribución conjunta de las estimaciones originales de los parámetros, entonces es fácil simular a partir de la solución OLS: simplemente genere realizaciones normales multivariadas de $\beta$ y convertirlos en los parámetros. Aquí se muestra una matriz de dispersión de 2000 realizaciones de este tipo. La fuerte curvilineidad muestra por qué el método Delta puede dar malos resultados.

Answer 2

Para ello existe un método estándar llamado método delta. Usted forma la inversa del hessiano de la log-verosimilitud wrt sus cuatro parámetros. Hay un parámetro extra para la varianza de de los residuos, pero no juega un papel en estos cálculos. A continuación, se calcula la respuesta predicha para los valores deseados de la variable independiente y se calcula su gradiente (la derivada con respecto a estos cuatro parámetros). estos cuatro parámetros. Llamamos a la inversa del hessiano $I$ y el vector gradiente $g$ . Se forma el producto vectorial matricial $$-g^tIg$$ Esto le da la varianza estimada para esa variable dependiente. Saque la raíz cuadrada para obtener la desviación estándar estimada. límites de confianza son el valor predicho +- dos desviaciones estándar. Para el caso especial de una regresión no lineal, puede corregir los grados de libertad. Usted tiene 10 observaciones y 4 parámetros, por lo que puede aumentar la estimación de la varianza en el modelo multiplicando por 10/6. Varios paquetes de software lo harán por usted. Escribí tu modelo en AD Model en AD Model Builder y lo ajusté y calculé las varianzas (sin modificar). Serán ligeramente diferentes de las suyas porque tuve que adivinar un poco los valores.

                    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

Esto puede hacerse para cualquier variable dependiente en el Constructor de Modelos AD. Se declara una variable en el lugar apropiado del código, así

   sdreport_number dep

y escribe el código para evaluar la variable dependiente así

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

Obsérvese que esto se evalúa para un valor de la variable independiente 2 veces el mayor observado en el ajuste del modelo. Si se ajusta el modelo, se obtiene la desviación estándar de esta variable dependiente

19   dep          7.2535e+00 1.0980e-01

He modificado el programa para incluir el código para calcular la límites de confianza para la función entalpía-volumen El archivo de código (TPL) tiene el siguiente aspecto

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)));

Luego volví a ajustar el modelo para obtener los desvíos estándar para las estimaciones de 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

Se han calculado para los valores de V observados, pero podrían calcularse fácilmente para cualquier valor de V.

Se ha señalado que en realidad se trata de un modelo lineal para el que existe un sencillo código R para realizar la estimación de los parámetros mediante OLS. Esto es muy atractivo, especialmente para los usuarios ingenuos. Sin embargo, desde el trabajo de Huber, hace más de treinta años, sabemos, o deberíamos saber, que casi siempre hay que sustituir OLS por una alternativa moderadamente robusta. La razón por la que esto no se hace de forma rutinaria creo que es que los métodos robustos son inherentemente no lineales. Desde este punto de vista, los métodos OLS simples y atractivos en R son más una trampa que una característica. Una ventaja del enfoque del Constructor de Modelos AD es su soporte incorporado para la modelización no lineal. Para cambiar el código de mínimos cuadrados a una mezcla normal robusta, sólo hay que cambiar una línea del código línea del código. La línea

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

se cambia a

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))));

La cantidad de sobredispersión en los modelos se mide por el parámetro a. Si a es igual a 1,0, la varianza es la misma que la del modelo normal. Si la varianza está inflada por los valores atípicos, se espera que a sea inferior a 1,0. Para estos datos, la estimación de a es de aproximadamente 0,23, por lo que la varianza es aproximadamente 1/4 de la varianza del modelo normal. La interpretación es que los valores atípicos han aumentado la estimación de la varianza en un factor de aproximadamente 4. El efecto de esto es aumentar el tamaño de los límites de confianza para los parámetros del modelo OLS. Esto representa una pérdida de eficiencia. Para el modelo de mezcla normal, las desviaciones estándar estimadas para la función entalpía-volumen son

 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

Se observa que hay pequeños cambios en las estimaciones puntuales, mientras que los límites de confianza se han reducido a cerca del 60% de los producidos por OLS.

El punto principal que quiero hacer es que todos los cálculos modificados ocurren automáticamente una vez que uno cambia la única línea de código en el archivo TPL.

Answer 3

La validación cruzada es una forma sencilla de estimación fiabilidad de su curva: https://en.wikipedia.org/wiki/Cross-validation_(estadísticas)

La propagación de la incertidumbre con diferenciales parciales es genial si realmente sabes $\Delta E_{0},\Delta V_{0},\Delta B_{0}$ y $\Delta B'$ . Sin embargo, el programa que está utilizando sólo da errores de ajuste (?). Esos serán demasiado optimistas (irrealmente pequeños).

Puede calcular el error de validación de una vez dejando uno de sus puntos fuera del ajuste y utilizando la curva ajustada para predecir el valor del punto que se dejó fuera. Repita esta operación para todos los puntos, de modo que cada uno de ellos se deje fuera una vez. A continuación, calcule el error de validación de su curva final (curva ajustada con todos los puntos) como media de los errores de predicción.

Esto sólo le dirá lo sensible que es su modelo para cualquier punto de datos nuevo. Por ejemplo, no le dirá lo inexacto que es su modelo de energía. Sin embargo, esto será mucho más realista estimación de error mero error de ajuste.

Además, si lo desea, puede trazar los errores de predicción en función del volumen.