Wie aussagekräftig oder nützlich sind Parameterfehler, wenn ein ungewichtetes LinearModelFit oder NonlinearModelFit ausgeführt wird?

Aug 16 2020

Dies mag eine Frage sein, die kurz davor steht, mehr zum Bereich der Statistik und der Cross-Validated SE zu gehören , aber ich bin auch speziell an Mathematica-Anpassungsroutinen interessiert.

Wenn ich ein Modell mit einem NonlinearModelFitoder mehreren Daten an einige Daten anpassen möchte, sind normalerweise LinearModelFiteinige Fehler mit meinem verbunden$y$-Daten, mit denen ich die Passungen gewichte. Diese Fehler können einfach der Standardfehler sein, der durch wiederholte Messungen erhalten wurde, oder ich weiß etwas über die physikalischen Prozesse und kann Gewichte zuweisen.

Zum Beispiel Weights->1/YDataErrors^2und ich setze meinen Varianzschätzer immer auf VarianceEstimatorFunction -> (1 &). Ich kann dann meine Parameterfehler aus der Kovarianzmatrix oder einfach mit abrufen MyFit["ParameterErrors"].

In einigen Fällen treten jedoch möglicherweise keine Fehler für die Daten auf, die angepasst werden sollen, was bedeutet, dass keine Gewichte wie oben beschrieben angegeben werden können. Meine Frage ist dann, wie zuverlässig - oder was noch wichtiger ist - wie physikalisch / statistisch bedeutsam sind die Parameterfehler für eine ungewichtete Anpassung in Mathematica?

Antworten

4 JimB Aug 16 2020 at 19:20

Wenn man beispielsweise zwei Fehlerquellen hat, beispielsweise einen Messfehler und einen Fehlanpassungsfehler, kann die Verwendung der auf den Messfehlern basierenden Gewichte zu einer starken Unterschätzung der Standardfehler führen. Betrachten Sie das folgende Modell:

$$y=a+b x +\gamma + \epsilon$$

wo $y$ ist die gemessene Antwort, $x$ ist der Prädiktor, $a$ und $b$ sind zu schätzende Konstanten, $\gamma$ ist der wiederholte Messfehler mit $\gamma \sim N(0,\sigma_{ME})$, und $\epsilon$ ist der Fehlanpassungsfehler mit $\epsilon \sim N(0,\sigma)$ und alle Fehler werden als unabhängig angenommen.

Stellen Sie zuerst einige spezifische Parameter ein:

(* Measurement error standard deviation *)
σME = 10;

(* Lack-of-fit error standard deviation *)
σ = 20;

(* Regression coefficients *)
a = 1;
b = 1;

Generieren und zeichnen Sie einige Daten:

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]

Betrachten Sie nun zwei verschiedene lineare Modellanpassungen, wo lm1Sie dies vorschlagen und lm2was ich vorschlage:

lm1 = LinearModelFit[data, z, z, Weights -> 1/ConstantArray[σME^2, n],
   VarianceEstimatorFunction -> (1 &)];
lm2 = LinearModelFit[data, z, z];
lm1["ParameterTable"]

lm2["ParameterTable"]

Die Schätzungen der Parameter sind identisch, aber die Standardfehler für lm1sind weniger als halb so groß wie für lm2. Was ist richtig?

Die "wahre" Kovarianzmatrix der Schätzer der kleinsten Quadrate von aund bfür dieses Modell ist

$$\left(\sigma ^2+\sigma_{ME}^2\right) \left(X^T.X\right)^{-1}$$

wo $X$ist die Designmatrix. In Mathematica Code der Standardfehler für bIS

X = Transpose[{ConstantArray[1, n], Range[n]}]
Sqrt[(σME^2 + σ^2) Inverse[Transpose[X].X][[2, 2]]] // N
(* 0.0774635 *)

Das passt ziemlich gut dazu lm2.

Dies ist insofern ein leicht erfundenes Beispiel, als ich alle Messstandardfehler identisch habe, da die Regressionsfunktionen von Mathematica nur einen einzigen Fehlerterm zulassen. Wenn die Messstandardfehler identisch sind, ergibt sich ein äquivalentes Modell mit einem einzelnen Fehler.

Selbst wenn die Messstandardabweichungen erheblich variieren, bleibt das Problem der falschen Gewichtung bestehen, so dass sie nicht mit der Fehlerstruktur des Modells übereinstimmt.

Die Regressionsroutinen von Mathematica sind für Modelle mit mehr als einer Fehlerquelle noch nicht ausreichend. Ich wünschte, sie wären es.