Zuverlässigkeit einer angepassten Kurve?

Ich möchte die Unsicherheit oder Zuverlässigkeit einer angepassten Kurve abschätzen. Ich nenne absichtlich keine genaue mathematische Größe, nach der ich suche, da ich nicht weiß, was es ist.

Hier ist (Energie) die abhängige Variable (Antwort) und (Volumen) die unabhängige Variable. Ich möchte die Energie-Volumen-Kurve eines Materials finden. Also habe ich einige Berechnungen mit einem quantenchemischen Computerprogramm durchgeführt, um die Energie für einige Probenvolumina zu erhalten (grüne Kreise im Diagramm). $E$ $V$ $E(V)$

Dann habe ich diese Datenproben mit der Birch-Murnaghan-Funktion ausgestattet : was von vier Parametern abhängt: . Ich gehe auch davon aus, dass dies die richtige Anpassungsfunktion ist, sodass alle Fehler nur vom Rauschen der Samples herrühren. Im Folgenden wird die angepasste Funktion wird als Funktion der geschrieben werden .

E. (E. | V.) = {E.}_{0} + \frac{9 {V.}_{0} {B.}_{0}}{16} {{[{(\frac{{V.}_{0}}{V.})}^{\frac{2}{3}} - - 1]]}^{3} {B.}_{0}^{'} + {[{(\frac{{V.}_{0}}{V.})}^{\frac{2}{3}} - - 1]]}^{2} [6 - - 4 {(\frac{{V.}_{0}}{V.})}^{\frac{2}{3}}]]}},

$\mathbb{E}(E|V) = E_0 + \frac{9V_0B_0}{16} \left\{ \left[\left(\frac{V_0}{V}\right)^\frac{2}{3}-1\right]^3B_0^\prime + \left[\left(\frac{V_0}{V}\right)^\frac{2}{3}-1\right]^2 \left[6-4\left(\frac{V_0}{V}\right)^\frac{2}{3}\right]\right\}\;,$

E_{0}, V_{0}, B_{0}, B_{0}^{'}

$E_0, V_0, B_0, B_0'$

(\hat{E})

$(\hat{E})$

V

$V$

Hier sehen Sie das Ergebnis (Anpassung mit einem Algorithmus der kleinsten Quadrate). Das y-Achsen - Variable ist und der X-Achsen - Variable . Die blaue Linie ist die Anpassung und die grünen Kreise sind die Stichprobenpunkte. $E$ $V$

Ich brauche jetzt ein Maß für die Zuverlässigkeit (am besten in Abhängigkeit des Volumens) diese angepaßten , weil ich es muß weitere Mengen berechnen , wie Übergangsdrücke oder Enthalpien. $\hat{E}(V)$

Meine Intutition sagt mir, dass die angepasste Kurve in der Mitte am zuverlässigsten ist, daher denke ich, dass die Unsicherheit (z. B. der Unsicherheitsbereich) gegen Ende der Probendaten zunehmen sollte, wie in dieser Skizze:

Was für eine Art von Maß suche ich jedoch und wie kann ich es berechnen?

Um genau zu sein, gibt es hier tatsächlich nur eine Fehlerquelle: Die berechneten Abtastwerte sind aufgrund von Rechengrenzen verrauscht. Wenn ich also einen dichten Satz von Datenproben berechnen würde, würden sie eine holprige Kurve bilden.

Meine Idee, die gewünschte Unsicherheitsschätzung zu finden, besteht darin, den folgenden "Fehler" basierend auf den Parametern zu berechnen, die Sie in der Schule lernen ( Ausbreitung der Unsicherheit ):

Dieundwerden von der Anpassungssoftware angegeben.

Δ E. (V.) = \sqrt{{(\frac{\partial E. (V.)}{\partial {E.}_{0}} Δ {E.}_{0})}^{2} + {(\frac{\partial E. (V.)}{\partial {V.}_{0}} Δ {V.}_{0})}^{2} + {(\frac{\partial E. (V.)}{\partial {B.}_{0}} Δ {B.}_{0})}^{2} + {(\frac{\partial E. (V.)}{\partial {B.}_{0}^{'}} Δ {B.}_{0}^{'})}^{2}}

$\Delta E(V) = \sqrt{ \left(\frac{\partial E(V)}{\partial E_0} \Delta E_0\right)^2 + \left(\frac{\partial E(V)}{\partial V_0} \Delta V_0\right)^2 + \left(\frac{\partial E(V)}{\partial B_0} \Delta B_0\right)^2 + \left(\frac{\partial E(V)}{\partial B_0'} \Delta B_0'\right)^2}$

Δ E_{0}, Δ V_{0}, Δ B_{0}

$\Delta E_0, \Delta V_0, \Delta B_0$

Δ B_{0}^{'}

$\Delta B_0'$

Ist das ein akzeptabler Ansatz oder mache ich es falsch?

PS: Ich weiß, dass ich auch nur die Quadrate der Residuen zwischen meinen Datenproben und der Kurve zusammenfassen könnte, um eine Art Standardfehler zu erhalten, aber dies ist nicht volumenabhängig.

error fitting reliability uncertainty Thymian
quelle

Keiner Ihrer Parameter ist ein Exponent, was gut ist. Welche NLS-Software haben Sie verwendet? Die meisten geben eine Schätzung für die parametrische Unsicherheit zurück (was völlig unrealistisch sein kann, wenn Ihre Parameter Exponenten sind, dies ist jedoch nicht der Fall).

DeltaIV

Es gibt kein A auf der rechten Seite Ihrer Gleichung, aber es erscheint in Ihrer Darstellung. Wenn Sie "vier Parameter" sagen, meinen Sie Parameter im statistischen Sinne (in welchem Fall, wo sind Ihre IVs) oder meinen Sie Variablen (in welchem Fall, wo sind Ihre Parameter)? Bitte klären Sie die Rollen der Symbole - was wird gemessen und was sind Unbekannte?

Glen_b -State Monica

Ich denke das V ist A ^ 3. Das habe ich benutzt und meine Handlung sah identisch mit seiner aus.

Dave Fournier

@Glen_b Ich habe gerade angenommen, dass die Y-Achse in der Birch-Murnaghan-Funktion E ist, während die x-Achse V ist. Die vier Parameter sind die vier Parameter in der Birch-Murnaghan-Funktion. Wenn Sie davon ausgehen, dass Sie etwas bekommen, das so aussieht, wie er es hat.

Dave Fournier

Ah, warte, ich verstehe es endlich.

ist kein Erwartungsoperator (wie ich auf der LHS einer Gleichung ohne Fehlerterm auf der RHS erwarten würde),

ist die Antwortvariable, die als Funktion in der Form

. GROSSER HINWEIS für alle: Zeigen Sie einem Statistiker keine Gleichung mit

links von einer Regressionsgleichung, ohne genau zu definieren, was Sie meinen, da er wahrscheinlich davon ausgeht, dass es sich um eine Erwartung handelt.

E ()

$E()$

E

$E$

y (x)

$y(x)$

E ()

$E()$

Glen_b -Reinstate Monica

Antworten:

Dies ist ein gewöhnliches Problem der kleinsten Quadrate!

Definieren

x = V^{- 2 /. 3}, w = {V.}_{0}^{1 /. 3},

$x = V^{-2/3}, \ w = V_0^{1/3},$

Das Modell kann neu geschrieben werden

E. (E. | V.) = β_{0} + β_{1} x + β_{2} x^{2} + β_{3} x^{3}

$\mathbb{E}(E|V) = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3$

$\beta=(\beta_i)^\prime$

16 β = (\begin{matrix} 16 {E.}_{0} + 54 {B.}_{0} w^{3} - - 9 {B.}_{0} {B.}_{0}^{'} w^{3} \\ - - 144 {B.}_{0} w^{5} + 27 {B.}_{0} {B.}_{0}^{'} w^{5} \\ 126 {B.}_{0} w^{7} - - 27 {B.}_{0} {B.}_{0}^{'} w^{7} \\ - - 36 {B.}_{0} w^{9} + 9 {B.}_{0} {B.}_{0}^{'} w^{9} \end{matrix}) .

$16 \beta = \pmatrix{16 E_0 + 54 B_0 w^3 - 9 B_0 B_0^\prime w^3\\ -144 B_0 w^5 + 27 B_0 B_0^\prime w^5\\ 126 B_0 w^7 - 27 B_0 B_0^\prime w^7\\ -36 B_0 w^9 + 9 B_0 B_0^\prime w^9}.$

$B_0, B_0^\prime$ $w$ $B_0, B_0^\prime, w$ $E_0$ $\beta$

$(E_0, B_0, B_0^\prime, V_0)$ $E$

$\hat\beta$ R

#
# The data.
#
X <- data.frame(V=c(41, 43, 46, 48, 51, 53, 55.5, 58, 60, 62.5),
                E=c(-48.05, -48.5, -48.8, -49.03, -49.2, -49.3, -49.35, 
                    -49.34, -49.31, -49.27))
#
# OLS regression.
#
fit <- lm(E ~ I(V^(-2/3)) + I(V^(-4/3)) + I(V^(-6/3)), data=X)
summary(fit)
beta <- coef(fit)
#
# Prediction, including standard errors of prediction.
#
V0 <- seq(40, 65)
y <- predict(fit, se.fit=TRUE, newdata=data.frame(V=V0))
#
# Plot the data, the fit, and a three-SEP band.
#
plot(X$V, X$E, xlab="Volume", ylab="Energy", bty="n", xlim=c(40, 60))
polygon(c(V0, rev(V0)), c(y$fit + 3*y$se.fit, rev(y$fit - 3*y$se.fit)),
        border=NA, col="#f0f0f0")
curve(outer(x^(-2/3), 0:3, `^`) %*% beta, add=TRUE, col="Red", lwd=2)
points(X$V, X$E)

$\beta$

whuber
quelle

Zwar sind Algorithmen zum Anpassen linearer Modelle numerisch viel stabiler als für nichtlineare Modelle, doch gibt es keinen Unterschied in der Genauigkeit der Diagnose, solange der nichtlineare Anpassungsalgorithmus konvergiert. Ich habe nachgesehen und wir haben die gleiche Restsumme von Quadraten zu mindestens 4 Sig Feigen. Auch die von Ihnen gewählte lineare Parametrisierung ist sehr verwirrend, so dass laut t-Test keiner der Parameter signifikant ist. Alle von mir sind. Nicht wirklich eine große Sache, aber amüsant und könnte den jungen Spieler verwirren.

Dave Fournier

Ich denke auch, dass Sie die Frage des OP nicht beantwortet haben, da sie erklärte, sie wolle so etwas wie Vertrauensgrenzen für die Enthalpie-Volumen-Funktion

Dave Fournier

β

$\beta$

(E_{0}, \dots)

$(E_0,\ldots)$

({\hat{E}}_{0} \dots)

$(\hat E_0\ldots)$

whuber

Ihr Modell und meins sind unabhängig von der Parametrierung identisch. (Ich spreche über das OLS-Modell.) Wenn ein bestimmter Parameter linear in das Modell eingeht, führen die Standardabweichungen zu besseren Konfidenzgrenzen für diesen Parameter. Die über die Delta-Methode erhaltenen Standardabweichungen sind gleich, unabhängig davon, ob sie zur Parametrisierung des Modells verwendet oder als abhängige Variable gelöst werden. In diesem Fall ist die abhängige interessierende Variable die Enthalpie-Volumen-Funktion und ihre Delta-Methode std dev ist dieselbe, unabhängig davon, ob Sie Ihre oder meine Parametrisierung verwenden.

Dave Fournier

\hat{β}

$\hat\beta$

$I$ $g$

- - G^{t} ich G

$-g^tIg$ Dies gibt Ihnen die geschätzte Varianz für diese abhängige Variable. Nehmen Sie die Quadratwurzel, um die geschätzte Standardabweichung zu erhalten. dann sind die Konfidenzgrenzen der vorhergesagte Wert + - zwei Standardabweichungen. Dies ist Standard Likelihood Zeug. Für den Sonderfall einer nichtlinearen Regression können Sie die Freiheitsgrade korrigieren. Sie haben 10 Beobachtungen und 4 Parameter, sodass Sie die Schätzung der Varianz im Modell durch Multiplikation mit 10/6 erhöhen können. Mehrere Softwarepakete erledigen dies für Sie. Ich habe Ihr Modell in AD Model in AD Model Builder geschrieben und angepasst und die (unveränderten) Abweichungen berechnet. Sie werden sich geringfügig von Ihren unterscheiden, da ich die Werte etwas erraten musste.

                    estimate   std dev
10   pred_E      -4.8495e+01 7.5100e-03
11   pred_E      -4.8810e+01 7.9983e-03
12   pred_E      -4.9028e+01 7.5675e-03
13   pred_E      -4.9224e+01 6.4801e-03
14   pred_E      -4.9303e+01 6.8034e-03
15   pred_E      -4.9328e+01 7.1726e-03
16   pred_E      -4.9329e+01 7.0249e-03
17   pred_E      -4.9297e+01 7.1977e-03
18   pred_E      -4.9252e+01 1.1615e-02

Dies kann für jede abhängige Variable in AD Model Builder durchgeführt werden. Man deklariert eine Variable an der entsprechenden Stelle im Code wie folgt

   sdreport_number dep

und schreibt den Code, um die abhängige Variable wie folgt auszuwerten

dep=sqrt(V0-cube(Bp0)/(1+2*max(V)));

Beachten Sie, dass dies für einen Wert der unabhängigen Variablen ausgewertet wird, der doppelt so groß ist wie der größte Wert, der in der Modellanpassung beobachtet wurde. Passen Sie das Modell an und Sie erhalten die Standardabweichung für diese abhängige Variable

19   dep          7.2535e+00 1.0980e-01

Ich habe das Programm so geändert, dass es Code zur Berechnung der Konfidenzgrenzen für die Enthalpie-Volumen-Funktion enthält. Die TPL-Datei (Code) sieht aus

DATA_SECTION
 init_int nobs
 init_matrix data(1,nobs,1,2)
 vector E
 vector V
 number Vmean
LOC_CALCS
 E=column(data,2);
 V=column(data,1);
 Vmean=mean(V);

PARAMETER_SECTION
 init_number E0
 init_number log_V0_coff(2)
 init_number log_B0(3)
 init_number log_Bp0(3)
 init_bounded_number a(.9,1.1)
 sdreport_number V0
 sdreport_number B0
 sdreport_number Bp0
 sdreport_vector pred_E(1,nobs)
 sdreport_vector P(1,nobs)
 sdreport_vector H(1,nobs)
 sdreport_number dep
 objective_function_value f
PROCEDURE_SECTION
  V0=exp(log_V0_coff)*Vmean;
  B0=exp(log_B0);
  Bp0=exp(log_Bp0);
  if (current_phase()<4)
  f+=square(log_V0_coff) +square(log_B0);

  dvar_vector sv=pow(V0/V,0.66666667);
  pred_E=E0 + 9*V0*B0*(cube(sv-1.0)*Bp0
    + elem_prod(square(sv-1.0),(6-4*sv)));

  dvar_vector r2=square(E-pred_E);
  dvariable vhat=sum(r2)/nobs;
  dvariable v=a*vhat;
  f=0.5*nobs*log(v)+sum(r2)/(2.0*v);

  // code to calculate the  enthalpy-volume function
  double delta=1.e-4;
  dvar_vector svp=pow(V0/(V+delta),0.66666667);
  dvar_vector svm=pow(V0/(V-delta),0.66666667);
  P = -((9*V0*B0*(cube(svp-1.0)*Bp0
      + elem_prod(square(svp-1.0),(6-4*svp))))
      -(9*V0*B0*(cube(svm-1.0)*Bp0
      + elem_prod(square(svm-1.0),(6-4*svm)))))/(2.0*delta);
  H=E+elem_prod(P,V);

dep=sqrt(V0-cube(Bp0)/(1+2*max(V)));

Dann habe ich das Modell umgerüstet, um die Standardentwickler für die Schätzungen von H zu erhalten.

29   H           -3.9550e+01 5.9163e-01
30   H           -4.1554e+01 2.8707e-01
31   H           -4.3844e+01 1.2333e-01
32   H           -4.5212e+01 1.5011e-01
33   H           -4.6859e+01 1.5434e-01
34   H           -4.7813e+01 1.2679e-01
35   H           -4.8808e+01 1.1036e-01
36   H           -4.9626e+01 1.8374e-01
37   H           -5.0186e+01 2.8421e-01
38   H           -5.0806e+01 4.3179e-01

Diese werden für Ihre beobachteten V-Werte berechnet, können jedoch leicht für jeden Wert von V berechnet werden.

Es wurde darauf hingewiesen, dass dies tatsächlich ein lineares Modell ist, für das es einen einfachen R-Code gibt, um die Parameterschätzung über OLS durchzuführen. Dies ist besonders für naive Benutzer sehr ansprechend. Seit der Arbeit von Huber vor über dreißig Jahren wissen wir jedoch oder sollten wissen, dass man OLS wahrscheinlich fast immer durch eine mäßig robuste Alternative ersetzen sollte. Der Grund, warum dies meiner Meinung nach nicht routinemäßig gemacht wird, ist, dass robuste Methoden von Natur aus nichtlinear sind. Unter diesem Gesichtspunkt sind die einfachen ansprechenden OLS-Methoden in R eher eine Falle als ein Merkmal. Ein Vorteil des AD Model Builder-Ansatzes ist die integrierte Unterstützung für nichtlineare Modellierung. Um den Code der kleinsten Quadrate in eine robuste normale Mischung zu ändern, muss nur eine Codezeile geändert werden. Die Linie

    f=0.5*nobs*log(v)+sum(r2)/(2.0*v);

wird geändert in

f=0.5*nobs*log(v)
  -sum(log(0.95*exp(-0.5*r2/v) + 0.05/3.0*exp(-0.5*r2/(9.0*v))));

Das Ausmaß der Überdispersion in den Modellen wird durch den Parameter a gemessen. Wenn a gleich 1,0 ist, ist die Varianz dieselbe wie für das normale Modell. Wenn die Varianz durch Ausreißer aufgeblasen wird, erwarten wir, dass a kleiner als 1,0 ist. Für diese Daten beträgt die Schätzung von a ungefähr 0,23, so dass die Varianz ungefähr 1/4 der Varianz für das normale Modell beträgt. Die Interpretation ist, dass Ausreißer die Varianzschätzung um einen Faktor von ungefähr 4 erhöht haben. Dies hat zur Folge, dass die Konfidenzgrenzen für Parameter für das OLS-Modell vergrößert werden. Dies bedeutet einen Effizienzverlust. Für das normale Mischungsmodell betragen die geschätzten Standardabweichungen für die Enthalpievolumenfunktion

 29   H           -3.9777e+01 3.3845e-01
 30   H           -4.1566e+01 1.6179e-01
 31   H           -4.3688e+01 7.6799e-02
 32   H           -4.5018e+01 9.4855e-02
 33   H           -4.6684e+01 9.5829e-02
 34   H           -4.7688e+01 7.7409e-02
 35   H           -4.8772e+01 6.2781e-02
 36   H           -4.9702e+01 1.0411e-01
 37   H           -5.0362e+01 1.6380e-01
 38   H           -5.1114e+01 2.5164e-01

Man sieht, dass sich die Punktschätzungen geringfügig ändern, während die Konfidenzgrenzen auf etwa 60% der von OLS erstellten Grenzwerte gesenkt wurden.

Der wichtigste Punkt, den ich ansprechen möchte, ist, dass alle geänderten Berechnungen automatisch erfolgen, sobald die eine Codezeile in der TPL-Datei geändert wird.

Dave Fournier
quelle

I

$I$

E (E ∣ V)

$\mathbb{E}(E\mid V)$

E (E ∣ V)

$\mathbb{E}(E\mid V)$

E (H ∣ V)

$\mathbb{E}(H \mid V)$

@jwimberley, Sie sagen im Grunde, dass Dave Fourier die Formel für das Konfidenzintervall des (bedingten) Mittelwerts angegeben hat, während Thymian möglicherweise am Vorhersageintervall für eine neue Beobachtung interessiert ist. Letzteres ist für OLS einfach zu berechnen. Wie berechnen Sie es in diesem Fall?

DeltaIV

E = f (V) + ϵ

$E = f(V) + \epsilon$

E - \hat{E}

$E -\hat E$

ϵ

$\epsilon$

V

$V$

ϵ

$\epsilon$

ϵ

$\epsilon$

Jwimberley

@jwimberley Ich habe nur die Konfidenzgrenzen für die vorhergesagten Werte angezeigt, die den beobachteten V-Werten entsprechen, nur weil sie verfügbar waren. Ich habe meine Antwort bearbeitet, um zu zeigen, wie Konfidenzgrenzen für abhängige Variablen ermittelt werden.

Dave Fournier

Die Kreuzvalidierung ist eine einfache Methode, um die Zuverlässigkeit Ihrer Kurve abzuschätzen : https://en.wikipedia.org/wiki/Cross-validation_(statistics)

$\Delta E_{0},\Delta V_{0},\Delta B_{0}$ $\Delta B'$

Sie können den 1-fachen Validierungsfehler berechnen, indem Sie einen Ihrer Punkte von der Anpassung fernhalten und die angepasste Kurve verwenden, um den Wert des weggelassenen Punkts vorherzusagen. Wiederholen Sie dies für alle Punkte, so dass jeder einmal weggelassen wird. Berechnen Sie dann den Validierungsfehler Ihrer endgültigen Kurve (Kurve mit allen Punkten) als Durchschnitt der Vorhersagefehler.

Hier erfahren Sie nur, wie empfindlich Ihr Modell für einen neuen Datenpunkt ist. Beispielsweise wird Ihnen nicht mitgeteilt, wie ungenau Ihr Energiemodell ist. Dies ist jedoch eine viel realistischere Fehlerschätzung, lediglich ein Anpassungsfehler.

Wenn Sie möchten, können Sie auch Vorhersagefehler als Funktion des Volumens darstellen.

Jman
quelle