¿Cuán significativos o útiles son los errores de parámetros que se producen al realizar un LinearModelFit o NonlinearModelFit no ponderado?
Esta puede ser una pregunta que se tambalea al borde de pertenecer más al ámbito de las estadísticas y SE con validación cruzada , pero también estoy específicamente interesado en las rutinas de ajuste de Mathematica.
Por lo general, si quiero ajustar un modelo a algunos datos usando NonlinearModelFit
o LinearModelFit
, tendré algunos errores asociados con mi$y$-datos que utilizo para ponderar los ajustes. Estos errores pueden ser simplemente el error estándar adquirido a partir de mediciones repetidas, o puedo saber algo sobre los procesos físicos y puedo asignar pesos.
Por ejemplo Weights->1/YDataErrors^2
, siempre configuro mi estimador de varianza como VarianceEstimatorFunction -> (1 &)
. Luego puedo obtener mis errores de parámetro de la matriz de covarianza, o simplemente con MyFit["ParameterErrors"]
.
Sin embargo, en algunos casos, es posible que no haya errores para los datos que desea ajustar, lo que significa que no se pueden proporcionar pesos de la forma que describí anteriormente. Entonces, mi pregunta es, ¿cuán confiables, o más importante aún, cuán significativos son física / estadísticamente los errores de parámetro para un ajuste no ponderado en Mathematica?
Respuestas
Por ejemplo, si uno tiene dos fuentes de error, digamos un error de medición y un error de falta de ajuste, entonces el uso de pesos basados en los errores de medición puede resultar en subestimaciones grandes de los errores estándar. Considere el siguiente modelo:
$$y=a+b x +\gamma + \epsilon$$
dónde $y$ es la respuesta medida, $x$ es el predictor, $a$ y $b$ son constantes a estimar, $\gamma$ es el error de medición repetido con $\gamma \sim N(0,\sigma_{ME})$y $\epsilon$ es el error de falta de ajuste con $\epsilon \sim N(0,\sigma)$ y se supone que todos los errores son independientes.
Primero configure algunos parámetros específicos:
(* Measurement error standard deviation *)
σME = 10;
(* Lack-of-fit error standard deviation *)
σ = 20;
(* Regression coefficients *)
a = 1;
b = 1;
Genere y grafique algunos datos:
n = 100;
x = Range[n];
SeedRandom[12345];
measurementError = RandomVariate[NormalDistribution[0, σME], n];
lackOfFitError = RandomVariate[NormalDistribution[0, σ], n];
y = a + b x + measurementError + lackOfFitError;
data = Transpose[{x, y}];
data2 = {#[[1]], Around[#[[2]], σME]} & /@ data;
ListPlot[data2]

Ahora considere dos ajustes de modelo lineal diferentes donde lm1
es lo que sugiere y lm2
es lo que sugiero:
lm1 = LinearModelFit[data, z, z, Weights -> 1/ConstantArray[σME^2, n],
VarianceEstimatorFunction -> (1 &)];
lm2 = LinearModelFit[data, z, z];
lm1["ParameterTable"]

lm2["ParameterTable"]

Las estimaciones de los parámetros son idénticas pero los errores estándar de lm1
son menos de la mitad del tamaño que los de lm2
. Cual es la correcta?
La matriz de covarianza "verdadera" de los estimadores de mínimos cuadrados de a
y b
para este modelo es
$$\left(\sigma ^2+\sigma_{ME}^2\right) \left(X^T.X\right)^{-1}$$
dónde $X$es la matriz de diseño. En el código de Mathematica , el error estándar para b
es
X = Transpose[{ConstantArray[1, n], Range[n]}]
Sqrt[(σME^2 + σ^2) Inverse[Transpose[X].X][[2, 2]]] // N
(* 0.0774635 *)
Eso coincide bastante bien con lm2
.
Este es un ejemplo un poco elaborado en el sentido de que tengo todos los errores estándar de medición idénticos porque las funciones de regresión de Mathematica solo permiten un solo término de error. Y al tener los errores estándar de medición idénticos, eso da como resultado un modelo equivalente con un solo error.
Sin embargo, incluso cuando las desviaciones estándar de la medición varían considerablemente, sigue existiendo el problema de la ponderación incorrecta, de modo que no coincida con la estructura de error del modelo.
Las rutinas de regresión de Mathematica aún no son adecuadas para modelos con más de una fuente de error. Ojalá lo fueran.