Solución:
Aquí hay un método para hacer una regresión ortogonal ponderada de una línea recta, basado en las fórmulas de Krystek / Anton y York:
ortlinfit[data_?MatrixQ, errs_?MatrixQ] :=
Module[n = Length[data], c, ct, dk, dm, k, m, p, s, st, ul, vl, w, wt, xm, ym,
(* yes, I know I could have used FindFit[] for this... *)
ct, st, k = Flatten[MapAt[Normalize[1, #] &,
NArgMin[Norm[Function[x, y, y - [FormalM] x - [FormalK]] @@@ data],
[FormalM], [FormalK]], 1]];
(* find orthogonal regression coefficients *)
c, s, p = FindArgMin[
Total[(data.-[FormalS], [FormalC] -
[FormalP])^2/((errs^2).[FormalS]^2, [FormalC]^2)],
[FormalC]^2 + [FormalS]^2 == 1,
[FormalC], ct, [FormalS], st, [FormalP], k/ct];
(* slope and intercept *)
m, k = s, p/c;
wt = 1/errs^2; w = (Times @@@ wt)/(wt.1, m^2);
xm, ym = w.data/Total[w];
ul = data[[All, 1]] - xm; vl = data[[All, 2]] - ym;
(* uncertainties in slope and intercept *)
dm = w.(m ul - vl)^2/w.ul^2/(n - 2);
dk = dm (w.data[[All, 1]]^2/Total[w]);
Function[[FormalX], Evaluate[m, k.[FormalX], 1]], Sqrt[dm, dk]] /;
Dimensions[data] === Dimensions[errs]
ortlinfit[] espera data para contener los $ (x_j, y_j) $ pares, y errs para contener las incertidumbres correspondientes $ ( rm dx _j, rm dy _j) $. La rutina devuelve la línea de mejor ajuste como una función pura, así como las incertidumbres en la pendiente y la intersección ($ sigma_m $ y $ sigma_k $).
Como ejemplo, aquí hay algunos datos utilizados en el artículo de York:
data = 0, 5.9, 0.9, 5.4, 1.8, 4.4, 2.6, 4.6, 3.3, 3.5, 4.4, 3.7,
5.2, 2.8, 6.1, 2.8, 6.5, 2.4, 7.4, 1.5;
errs = 1000., 1., 1000., 1.8, 500., 4., 800., 8., 200., 20.,
80., 20., 60., 70., 20., 70., 1.8, 100., 1, 500. // Sqrt[1/#] &;
lin, sm, sk = ortlinfit[data, errs]
Function[x, 5.47991 - 0.480533 x], 0.0710065, 0.361871
Ahora, veamos los datos, las elipses de error asociadas (construidas a partir de las incertidumbres), la línea de mejor ajuste y las “líneas delimitadoras” $ (m- sigma_m) x + (k- sigma_k) $ y $ (m + sigma_m) x + (k + sigma_k) $:
Show[
Graphics[AbsolutePointSize[4], Point[data], MapThread[Circle, data, errs],
Frame -> True],
Plot[lin[x], lin[x] - sm x - sk, lin[x] + sm x + sk, x, -1, 9,
PlotStyle -> Directive[Thick, Red],
Directive[Dashed, Gray], Directive[Dashed, Gray]]
]
Resumen
Utilizar el Weights opción a LinearModelFit, estableciendo los pesos para que sean inversamente proporcionales a las varianzas de los términos de error.
Teoría
Este es un problema estándar: cuando los errores en los valores individuales $ y $ se expresan de una manera que se puede relacionar con su variaciones, luego use mínimos cuadrados ponderados con las varianzas recíprocas como pesos. (Busque en nuestro sitio hermano Validado cruzado para obtener más información sobre esto, incluidas referencias y generalizaciones).
Creando datos realistas para estudiar
Para ilustrar, suponga que los datos se dan como vectores $ x $ y $ y $ con los “errores” expresados como desviaciones estándar de $ y $ o como errores estándar de estimación de $ y $, o cualquier otra cantidad que pueda interpretarse como un múltiplo constante fijo de las desviaciones estándar de $ y $. Específicamente, el modelo aplicable para estos datos es
$$ y_i = beta_0 + beta_1 x_i + varepsilon_i $$
donde $ beta_0 $ (la intersección) y $ beta_1 $ (la pendiente) son constantes a estimar, $ varepsilon_i $ son desviaciones aleatorias independientes con media cero y $ text Var ( varepsilon_i) = sigma_i ^ 2 $ para algunos dado cantidades $ sigma_i $ se supone que se conocen con precisión. (El caso en el que todos los $ sigma_i $ son iguales a una constante desconocida común es la configuración de regresión lineal habitual).
En Mathematica podemos simular dichos datos con generación de números aleatorios. Envuelva esto en una función cuyos argumentos son la cantidad de datos y la pendiente y la intersección. Haré que los tamaños de los errores varíen aleatoriamente, pero generalmente aumentarán a medida que aumente $ x $.
simulate[n_Integer, intercept_: 0, slope_: 0] /; n >= 1 :=
Module[x, y, errors, sim,
x = Range[n];
errors = RandomReal[GammaDistribution[n, #/(10 n)]] & /@ x;
y = RandomReal[NormalDistribution[intercept + slope x[[#]], errors[[#]]]] & / Range[n];
sim["x"] = x;
sim["y"] = y;
sim["errors"] = errors;
sim
]
Aquí hay un pequeño ejemplo de su uso.
SeedRandom[17];
x, y, errors = simulate[16, 50, -2/3][#] & /@ "x", "y", "errors";
ListPlot[y, y + errors, y - errors, Joined -> False, True, True,
PlotStyle -> PointSize[0.015], Thick, Thick,
AxesOrigin -> 0, Min[y - errors]]
Los puntos simulados están rodeados por bandas de error.
Estimación de mínimos cuadrados ponderados
Para ajustar estos datos, utilizar el Weights opción de LinearModelFit. Una vez más, preparémonos para un análisis posterior encapsulando el ajuste en una función. A modo de comparación, ajustemos los datos con y sin las ponderaciones.
trial[n_Integer: 1, intercept_: 0, slope_: 0] :=
Module[x, y, errors, t, fit, fit0,
x, y, errors = simulate[n, intercept, slope][#] & /@ "x", "y", "errors";
fit = LinearModelFit[y, t, t, Weights -> 1 / errors^2];
fit0 = LinearModelFit[y, t, t];
fit[#], fit0[#] & @ "BestFitParameters"
]
La salida es una lista cuyos elementos dan intersección, pendiente: el primer elemento es para el ajuste ponderado, el segundo para el no ponderado.
Comparación de Monte-Carlo de los métodos de mínimos cuadrados ponderados y ordinarios
Realicemos muchas pruebas independientes (digamos, $ 1000 $ de ellas para conjuntos de datos simulados de $ n = 32 $ puntos cada una) y comparemos sus resultados:
SeedRandom[17];
simulation = ParallelTable[trial[32, 20, -1/2], i, 1, 1000];
ranges = 18.5, 22, -0.65, -0.35;
TableForm[
Table[Histogram[simulation[[All, i, j]], ImageSize -> 250,
PlotRange -> ranges[[j]], Automatic, Axes -> True, False],
i, 1, 2, j, 1, 2],
TableHeadings -> "Weighted", "OLS", "Intercept", "Slope"
]
Debido a que especifiqué una intersección de $ 20 $ y una pendiente de $ -1 / 2 $, querremos usar estos valores como referencias. De hecho, los histogramas de la columna de la izquierda (“Intercepción”) muestran conjuntos de estimado intersecciones que rondan los $ 20 $ y los histogramas de la columna de la derecha (“Pendiente”) muestran conjuntos de pendientes que rondan los $ -0,50 = -1 / 2 $. Esto ilustra el hecho teórico de que las estimaciones en ambos casos son imparcial. Sin embargo, mirando más de cerca se extiende de los histogramas (lea los números en los ejes horizontales), vemos que los de la fila superior (“ponderados”) tener anchos más pequeños de sus contrapartes en la fila inferior (“MCO” o “Mínimos cuadrados ordinarios”). Esto es evidencia de que las estimaciones ponderadas tienden, en general, a ser mejor luego los no ponderados, porque tienden a desviarse menos del true valores paramétricos.
Cuando los datos subyacentes realmente se ajustan al modelo hipotetizado – existe una relación lineal entre los $ x $ y los $ y $, con errores en los $ y $ que tienen desviaciones estándar conocidas pero diferentes – entonces entre todas las estimaciones lineales insesgadas de la pendiente y la intersección, mínimos cuadrados ponderados que utilizan varianzas recíprocas para las ponderaciones son mejores en el sentido que se acaba de ilustrar.
Obtención del error de estimación de la pendiente
Ahora, para responder a la pregunta: nos gustaría evaluar el error de estimación en la pendiente. Esto se puede obtener de la fit objeto de muchas formas: consulte la página de ayuda para obtener más detalles. Aquí hay una tabla muy bien formateada:
fit = LinearModelFit[y, t, t, Weights -> 1 / errors^2];
fit0 = LinearModelFit[y, t, t];
TableForm[fit[#], fit0[#] & @ "ParameterConfidenceIntervalTable",
TableHeadings -> , "Weighted", "OLS"]
En este caso, para este conjunto particular de datos simulados (como se creó anteriormente), el método ponderado informa un error estándar mucho menor para el interceptar que el método OLS (porque los errores cercanos a $ x = 0 $ son bajos de acuerdo con la información en errors) pero la estimación ponderada de la Pendiente tiene solo un error estándar ligeramente menor que la estimación de la pendiente por MCO.
Comentarios
Los errores tanto en $ x $ como en $ y $ pueden manejarse utilizando, por ejemplo, métodos de máxima verosimilitud. Sin embargo, esto implica una maquinaria mucho más matemática, estadística y computacional y requiere una evaluación cuidadosa de la naturaleza de esos errores (por ejemplo, si los errores $ x $ y los errores $ y $ son independientes). Un resultado general en la literatura estadística es que cuando los errores $ x $ son típicamente más pequeños que los errores $ y $, pero independientes de ellos, generalmente es seguro ignorar los errores $ x $. Para obtener más información sobre todo esto, los buenos términos de búsqueda incluyen “regresión de errores en las variables”, “regresión de Deming” e incluso “análisis de componentes principales (PCA)”.
Hice esta implementación del método clásico (y fácil de entender) de York siguiendo este artículo de Cameron Reed.
f[x_List, y_List, wxi_List, wyi_List] :=
Module[n = [email protected], wi, ui, vi, wmean, d, g, a, b, set, least,
wi[i_, m_] := wxi[[i]] wyi[[i]]/(m ^2 wyi[[i]] + wxi[[i]]);
ui[i_, m_] := x[[i]] - wmean[x, m];
vi[i_, m_] := y[[i]] - wmean[y, m];
wmean[q_List, m_] := Sum[wi[i, m] q[[i]], i, n]/Sum[wi[i, m], i, n];
d[m_] := Sum[wi[i, m]^2 ui[i, m]^2/wxi[[i]], i, n];
g[m_] :=- Sum[wi[i, m] ui[i, m] vi[i, m], i, n]/d[m];
a[m_] :=2 Sum[wi[i, m]^2 ui[i, m] vi[i, m]/wxi[[i]], i, n]/(3 d[m]);
b[m_] := (Sum[wi[i, m]^2 vi[i, m]^2/wxi[[i]], i, n] -
Sum[wi[i, m] ui[i, m]^2, i, n])/(3 d[m]);
set = [email protected] Reduce[m^3 - 3 a[m] m m + 3 b[m] m - g[m] == 0 &&
c == wmean[y, m] - m wmean[x, m], m, c,
Backsubstitution -> True];
least = Sum[wxi[[i]] (x[[i]] - (y[[i]] - c)/m)^2 +
wyi[[i]] (y[[i]] - (m x[[i]] + c))^2, i, [email protected]] /.
set[[[email protected][m /. set, _Real]]];
set[[[email protected][m /. set, _Real]]][[Position[least, [email protected]][[1]]]]];
Uso
f[Range[10], 3 Range[10] + RandomReal[.2], Array[# &, 10], Array[# &, 10]]
(*
-> m -> 3., c -> 0.110805
*)