Jak znaczące lub przydatne są błędy parametrów podczas wykonywania nieważonego LinearModelFit lub NonlinearModelFit?
Może to być pytanie, które balansuje na krawędzi przynależności bardziej do sfery statystyk i Cross-Validated SE , ale jestem również szczególnie zainteresowany procedurami dopasowania Mathematica.
Zwykle, jeśli chcę dopasować model do niektórych danych przy użyciu jednego z nich, NonlinearModelFit
lub LinearModelFit
wystąpią błędy związane z plikiem$y$-dane których używam do ważenia pasowań. Błędy te mogą być po prostu błędem standardowym uzyskanym z powtarzanych pomiarów lub mogę wiedzieć coś o procesach fizycznych i mogę przypisać wagi.
Na przykład, Weights->1/YDataErrors^2
a mój estymator wariancji zawsze ustawiam jako VarianceEstimatorFunction -> (1 &)
. Mogę wtedy uzyskać błędy parametrów z macierzy kowariancji lub po prostu z MyFit["ParameterErrors"]
.
Jednak w niektórych przypadkach może nie być żadnych błędów dla danych, które chcesz dopasować, co oznacza, że nie można podać wag w sposób, który opisałem powyżej. Moje pytanie brzmi zatem, na ile wiarygodne - lub co ważniejsze - jak fizycznie / statystycznie znaczące są błędy parametrów dla dopasowania nieważonego w Mathematica?
Odpowiedzi
Na przykład, jeśli mamy dwa źródła błędów, powiedzmy błąd pomiaru i błąd niedopasowania, to użycie wag opartych na błędach pomiaru może spowodować rażące niedoszacowanie błędów standardowych. Rozważ następujący model:
$$y=a+b x +\gamma + \epsilon$$
gdzie $y$ to zmierzona odpowiedź, $x$ jest predyktorem, $a$ i $b$ są stałymi do oszacowania, $\gamma$ to powtarzający się błąd pomiaru z $\gamma \sim N(0,\sigma_{ME})$, i $\epsilon$ to błąd braku dopasowania z $\epsilon \sim N(0,\sigma)$ i zakłada się, że wszystkie błędy są niezależne.
Najpierw ustaw określone parametry:
(* Measurement error standard deviation *)
σME = 10;
(* Lack-of-fit error standard deviation *)
σ = 20;
(* Regression coefficients *)
a = 1;
b = 1;
Wygeneruj i wykreśl niektóre dane:
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]

Rozważ teraz dwa różne modele liniowe, w których lm1
sugerujesz, a lm2
ja proponuję:
lm1 = LinearModelFit[data, z, z, Weights -> 1/ConstantArray[σME^2, n],
VarianceEstimatorFunction -> (1 &)];
lm2 = LinearModelFit[data, z, z];
lm1["ParameterTable"]

lm2["ParameterTable"]

Szacunki parametrów są identyczne, ale standardowe błędy lm1
są o mniej niż połowę mniejsze niż w przypadku lm2
. Który jest prawidłowy?
„Prawdziwa” macierz kowariancji estymatorów najmniejszych kwadratów tego modelu a
i b
dla tego modelu to
$$\left(\sigma ^2+\sigma_{ME}^2\right) \left(X^T.X\right)^{-1}$$
gdzie $X$to macierz projektowa. W kodzie Mathematica standardowy błąd dla b
to
X = Transpose[{ConstantArray[1, n], Range[n]}]
Sqrt[(σME^2 + σ^2) Inverse[Transpose[X].X][[2, 2]]] // N
(* 0.0774635 *)
To całkiem dobrze pasuje do lm2
.
Jest to nieco wymyślony przykład, w którym wszystkie standardowe błędy pomiaru są identyczne, ponieważ funkcje regresji Mathematica dopuszczają tylko jeden składnik błędu. A mając identyczne standardowe błędy pomiaru, daje to równoważny model z jednym błędem.
Jednak nawet wtedy, gdy odchylenia standardowe pomiaru różnią się znacznie, pozostaje kwestia niewłaściwego ważenia, takiego, że nie pasuje ono do struktury błędów modelu.
Procedury regresji Mathematica nie są jeszcze odpowiednie dla modeli z więcej niż jednym źródłem błędów. Chciałabym, żeby tak było.