Jak znaczące lub przydatne są błędy parametrów podczas wykonywania nieważonego LinearModelFit lub NonlinearModelFit?

Aug 16 2020

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, NonlinearModelFitlub LinearModelFitwystą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^2a 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

4 JimB Aug 16 2020 at 19:20

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 lm1sugerujesz, a lm2ja 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 lm1są 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 ai bdla 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 bto

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.