Fiabilité d'une courbe ajustée?

11

Je voudrais estimer l'incertitude ou la fiabilité d'une courbe ajustée. Je ne nomme pas intentionnellement une quantité mathématique précise que je recherche, car je ne sais pas ce que c'est.

Ici, (énergie) est la variable dépendante (réponse) et (volume) est la variable indépendante. Je voudrais trouver la courbe énergie-volume, , de certains matériaux. J'ai donc fait quelques calculs avec un programme informatique de chimie quantique pour obtenir l'énergie pour certains volumes d'échantillons (cercles verts dans l'intrigue).EVE(V)

J'ai ensuite ajusté ces échantillons de données avec la fonction Birch – Murnaghan : qui dépend de quatre paramètres: . Je suppose également que c'est la fonction d'ajustement correcte, donc toutes les erreurs viennent simplement du bruit des échantillons. Dans ce qui suit, la fonction ajustée sera écrit en fonction de .

E(E|V)=E0+9V0B016{[(V0V)231]3B0+[(V0V)231]2[64(V0V)23]},
E0,V0,B0,B0(E^)V

Ici vous pouvez voir le résultat (ajustement avec un algorithme des moindres carrés). La variable de l' axe des y est et la variable x est l' axe . La ligne bleue est l'ajustement et les cercles verts sont les points d'échantillonnage.EV

Ajustement bouleau – Murnaghan (bleu) de l'échantillon (vert)

J'ai maintenant besoin d'une certaine mesure de la fiabilité (au mieux en fonction du volume) de cette courbe ajustée, , car j'en ai besoin pour calculer d'autres quantités comme les pressions de transition ou les enthalpies.E^(V)

Mon intuition me dit que la courbe ajustée est la plus fiable au milieu, donc je suppose que l'incertitude (disons la plage d'incertitude) devrait augmenter vers la fin des données de l'échantillon, comme dans ce croquis: entrez la description de l'image ici

Cependant, quel est ce type de mesure que je recherche et comment puis-je le calculer?

Pour être précis, il n'y a en fait qu'une seule source d'erreur ici: les échantillons calculés sont bruyants en raison des limites de calcul. Donc, si je calculais un ensemble dense d'échantillons de données, ils formeraient une courbe cahoteuse.

Mon idée pour trouver l'estimation d'incertitude souhaitée est de calculer l '' 'erreur' 'suivante en fonction des paramètres au fur et à mesure que vous l'apprenez à l'école ( propagation de l'incertitude ):

ΔE(V)=(E(V)E0ΔE0)2+(E(V)V0ΔV0)2+(E(V)B0ΔB0)2+(E(V)B0ΔB0)2
Les et , sont fournis par le logiciel d'ajustement.ΔE0,ΔV0,ΔB0ΔB0

Est-ce une approche acceptable ou je me trompe?

PS: Je sais que je pourrais aussi simplement résumer les carrés des résidus entre mes échantillons de données et la courbe pour obtenir une sorte d '«erreur standard», mais cela ne dépend pas du volume.

thym
la source
aucun de vos paramètres n'est un exposant, ce qui est bien. Quel logiciel NLS avez-vous utilisé? La plupart renverront une estimation de l'incertitude paramétrique (qui peut être complètement irréaliste si vos paramètres sont des exposants, mais ce n'est pas votre cas).
DeltaIV
Il n'y a pas de A sur le côté droit de votre équation, mais il apparaît dans votre graphique. Quand vous dites «quatre paramètres», voulez-vous dire des paramètres au sens statistique (dans ce cas, où sont vos IV) ou voulez-vous dire des variables (dans quel cas, où sont vos paramètres)? Veuillez clarifier les rôles des symboles - qu'est-ce qui est mesuré et quelles sont les inconnues?
Glen_b -Reinstate Monica
1
Je pense que le V est A ^ 3. c'est ce que j'ai utilisé et mon intrigue était identique à la sienne.
dave fournier
@Glen_b Je viens de supposer que l'axe Y est E dans la fonction Birch – Murnaghan tandis que l'axe x est V. Les quatre paramètres sont les quatre paramètres dans la fonction Birch – Murnaghan. Si vous supposez que vous obtenez quelque chose qui ressemble à ce qu'il a.
dave fournier
Ah, attendez, je comprends enfin. n'est pas un opérateur d'attente (comme je m'attendrais à le voir sur le LHS d'une équation sans terme d'erreur sur le RHS), est la variable de réponse écrite comme une fonction sous la forme . GRAND CONSEIL à tout le monde: ne montrez pas une équation avec à gauche d'une équation de régression à un statisticien sans définir soigneusement ce que vous voulez dire, car ils supposeront probablement que c'est une attente. E()Ey(x)E()
Glen_b -Reinstate Monica

Réponses:

8

C'est un problème de moindres carrés ordinaire!

Définition

x=V2/3, w=V01/3,

le modèle peut être réécrit

E(E|V)=β0+β1x+β2x2+β3x3

où les coefficients sont algébriquement liés aux coefficients d'origine viaβ=(βi)

16β=(16E0+54B0w39B0B0w3144B0w5+27B0B0w5126B0w727B0B0w736B0w9+9B0B0w9).

C'est simple à résoudre algébriquement ou numériquement: choisissez la solution pour laquelle et sont positifs. La seule raison de le faire est d'obtenir des estimations de et et de vérifier qu'elles sont physiquement significatives. Toutes les analyses de l'ajustement peuvent être effectuées en termes de .B0,B0wB0,B0,wE0β

Cette approche est non seulement beaucoup plus simple que l'ajustement non linéaire, elle est également plus précise: la matrice de variance-covariance pour renvoyée par un ajustement non linéaire n'est qu'une approximation quadratique locale de la distribution d'échantillonnage de ces paramètres, alors que (pour les erreurs normalement réparties dans la mesure de , de toute façon) les résultats OLS ne sont pas des approximations.(E0,B0,B0,V0)E

Les intervalles de confiance, les intervalles de prédiction, etc. peuvent être obtenus de la manière habituelle sans avoir besoin de trouver ces valeurs: calculez-les en termes d'estimations et de leur matrice variance-covariance. (Même Excel pourrait le faire!) Voici un exemple, suivi du code (simple) qui l'a produit.β^R

Figure

#
# 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)

Si vous êtes intéressé par la distribution conjointe des estimations de paramètres d'origine, il est facile de simuler à partir de la solution OLS: générez simplement des réalisations normales multivariées de et convertissez-les en paramètres. Voici une matrice de nuage de points de 2000 réalisations. La forte curvilinéarité montre pourquoi la méthode Delta est susceptible de donner de mauvais résultats.β

Figure 2

whuber
la source
1
S'il est vrai que les algorithmes d'ajustement des modèles linéaires sont beaucoup plus stables numériquement que ceux des modèles non linéaires, il n'est pas vrai qu'il existe une différence dans la précision des diagnostics tant que l'algorithme d'ajustement non linéaire converge. J'ai vérifié et nous avons la même somme résiduelle de carrés pour au moins 4 figues sig. De plus, le paramétrage linéaire que vous avez choisi est très confus de sorte qu'aucun des paramètres n'est significatif selon le test t. Tous les miens le sont. Pas vraiment un gros problème, mais amusant et pourrait dérouter le jeune joueur.
dave fournier
De plus, je suppose que vous n'avez pas répondu à la question de l'OP puisqu'elle a déclaré qu'elle voulait quelque chose comme des limites de confiance pour la fonction enthalpie-volume
Dave Fournier
1
@Dave Ce sont des points réfléchis, merci. Dans un problème physique, on ne s'intéresse généralement pas à la signification: la théorie implique toutes les variables. On se préoccupe plutôt d' estimer les valeurs. Bien que les deux approches devraient atteindre la même perte minimale (somme des carrés des résidus), OLS produit des distributions correctes pour la variance d'échantillonnage de ses paramètres. L'approche non linéaire ne fonctionne pas. Il est exact d'appliquer la transformation de la distribution de à , mais l'utilisation de covariances de n'est qu'une approximation. β(E0,)(E^0)
whuber
Votre modèle et le mien sont identiques indépendamment du paramétrage. (Je parle du modèle OLS.) Il est vrai que si un paramètre particulier entre dans le modèle linéairement, les écarts-types produisent de meilleures limites de confiance pour ce paramètre. les écarts-types obtenus via la méthode delta seront les mêmes, qu'ils soient utilisés pour paramétrer le modèle ou résolus comme variable dépendante. Dans ce cas, la variable dépendante d'intérêt est la fonction enthalpie-volume et sa méthode delta std dev sera la même, que l'on utilise votre paramétrage ou le mien.
dave fournier
1
@Dave Je suis d'accord: les modèles sont les mêmes mais avec des paramètres différents. Les avantages de la formulation OLS s'étendent au fait que les erreurs normales dans la réponse se traduisent par une solution exacte pour la distribution des estimations des paramètres , qui est facilement transformée en une estimation de la distribution quadridimensionnelle complète des estimations des paramètres d'origine. Bien que cela puisse être fait en utilisant le paramétrage d'origine, cela nécessiterait plus de travail numérique. De plus, les avantages conceptuels de voir le modèle d'origine n'est que du SLO déguisé sont importants. β^
whuber
3

Il existe une approche standard pour cela appelée la méthode delta. Vous formez l'inverse de la Hesse de la log-vraisemblance par rapport à vos quatre paramètres. Il existe un paramètre supplémentaire pour la variance des résidus, mais il ne joue aucun rôle dans ces calculs. Ensuite, vous calculez la réponse prévue pour les valeurs souhaitées de la variable indépendante et calculez son gradient (la dérivée par rapport à ces quatre paramètres). Appelons l'inverse de la Hesse et du vecteur de gradient . Vous produit matriciel vectoriel Ig

gtIg
Cela vous donne la variance estimée pour cette variable dépendante. Prenez la racine carrée pour obtenir l'écart type estimé. alors les limites de confiance sont la valeur prédite + - deux écarts-types. C'est un truc de vraisemblance standard. pour le cas particulier d'une régression non linéaire, vous pouvez corriger les degrés de liberté. Vous disposez de 10 observations et de 4 paramètres pour augmenter l'estimation de la variance dans le modèle en multipliant par 10/6. Plusieurs progiciels le feront pour vous. J'ai rédigé votre modèle dans AD Model dans AD Model Builder, je l'ai ajusté et j'ai calculé les variances (non modifiées). Ils seront légèrement différents des vôtres car je devais deviner un peu les valeurs.
                    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

Cela peut être fait pour n'importe quelle variable dépendante dans AD Model Builder. On déclare une variable à l'endroit approprié dans le code comme ceci

   sdreport_number dep

et écrit le code évaluer la variable dépendante comme ceci

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

Notez que cela est évalué pour une valeur de la variable indépendante 2 fois la plus grande observée dans l'ajustement du modèle. Ajuster le modèle et on obtient l'écart type pour cette variable dépendante

19   dep          7.2535e+00 1.0980e-01

J'ai modifié le programme pour inclure le code de calcul des limites de confiance pour la fonction enthalpie-volume Le fichier de code (TPL) ressemble à

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)));

J'ai ensuite réaménagé le modèle pour obtenir les devs standard pour les estimations de H.

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

Ceux-ci sont calculés pour vos valeurs V observées, mais pourraient être facilement calculés pour n'importe quelle valeur de V.

Il a été souligné qu'il s'agit en fait d'un modèle linéaire pour lequel il existe un code R simple pour effectuer l'estimation des paramètres via OLS. C'est très attrayant en particulier pour les utilisateurs naïfs. Cependant, depuis les travaux de Huber il y a plus de trente ans, nous savons ou devons savoir que l'on devrait probablement presque toujours remplacer OLS par une alternative modérément robuste. La raison pour laquelle cela n'est pas systématiquement fait, je crois, c'est que les méthodes robustes sont intrinsèquement non linéaires. De ce point de vue, les méthodes OLS attrayantes simples dans R sont plus un piège, plutôt qu'une fonctionnalité. Un avantage de l'approche AD Model Builder est sa prise en charge intégrée pour la modélisation non linéaire. Pour changer le code des moindres carrés en un mélange normal robuste, une seule ligne du code doit être modifiée. La ligne

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

est changé en

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))));

La quantité de surdispersion dans les modèles est mesurée par le paramètre a. Si a est égal à 1,0, la variance est la même que pour le modèle normal. S'il y a inflation de la variance par les valeurs aberrantes, nous nous attendons à ce que a soit inférieur à 1,0. Pour ces données, l'estimation de a est d'environ 0,23, de sorte que la variance est d'environ 1/4 de la variance pour le modèle normal. L'interprétation est que les valeurs aberrantes ont augmenté l'estimation de la variance d'un facteur d'environ 4. Cela a pour effet d'augmenter la taille des limites de confiance pour les paramètres du modèle OLS. Cela représente une perte d'efficacité. Pour le modèle de mélange normal, les écarts-types estimés pour la fonction enthalpie-volume sont

 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

On constate qu'il y a de petits changements dans les estimations ponctuelles, tandis que les limites de confiance ont été réduites à environ 60% de celles produites par OLS.

Le point principal que je veux faire est que tous les calculs modifiés se produisent automatiquement une fois que l'on change la seule ligne de code dans le fichier TPL.

dave fournier
la source
2
Pour le bénéfice de @ thym, je voudrais noter que la «méthode delta» est essentiellement la même procédure que la «propagation de l'incertitude» qui lui est familière et qui est communément enseignée aux scientifiques - au moins, ils sont la même procédure lorsque celle-ci est effectuée correctement. Une mise en garde est que la formule publiée dans la question ignore les corrélations entre les valeurs estimées des paramètres. Ignorer les corrélations revient à ne considérer que les éléments diagonaux de dans la méthode delta. I
jwimberley
1
Aussi pour @thyme, la propagation des incertitudes / la méthode delta ne produit que l'incertitude dans . En plus de cela, tout biais / variance dû au bruit. Je suppose que vous faites des prédictions sur des échantillons physiques, dont l'énergie / l'enthalpie / d'autres quantités thermodynamiques n'ont pas le même bruit que dans votre logiciel de simulation, mais dans les deux cas, cela ajoute une source de variance supplémentaire en plus de la variance dans ou qui est due à l'incertitude de l'ajustement. E(EV)E(EV)E(HV)
jwimberley
1
@jwimberley, vous dites essentiellement que Dave Fourier a donné la formule de l'intervalle de confiance de la moyenne (conditionnelle), tandis que le thym peut être intéressé par l'intervalle de prédiction pour une nouvelle observation. Il est facile de calculer ce dernier pour OLS. Comment le calculez-vous dans ce cas?
DeltaIV
1
@DeltaIV Ce serait toujours facile dans ce cas - si le modèle des moindres carrés non linéaire est correct, alors les résidus de l'ajustement sont distribués comme . Ainsi, la variance supplémentaire dans l'intervalle de prédiction est la variance des résidus d'ajustement. Ceci est lié à l'idée dans le post-script de la réponse (qui n'est pas dépendant du car le modèle d'ajustement n'est pas hétéroscédastique). Cependant, plus important encore, dans l'ajustement provient de limites de calcul, tandis que dans le monde provient de fluctuations thermodynamiques, qui ne sont probablement pas comparables. E=f(V)+ϵEE^ϵVϵϵ
jwimberley
1
@jwimberley I n'a montré que les limites de confiance pour les valeurs prévues correspondant aux valeurs V observées simplement parce qu'elles étaient disponibles. J'ai modifié ma réponse pour montrer comment obtenir des limites de confiance pour n'importe quelle variable dépendante.
dave fournier
0

La validation croisée est un moyen simple d' estimer la fiabilité de votre courbe: https://en.wikipedia.org/wiki/Cross-validation_(statistics)

La propagation de l'incertitude avec des différentiels partiels est excellente si vous connaissez vraiment et . Cependant, le programme que vous utilisez ne donne que des erreurs d'ajustement (?). Celles-ci seront trop optimistes (irréalistes).ΔE0,ΔV0,ΔB0ΔB

Vous pouvez calculer une erreur de validation 1 fois en laissant un de vos points loin de l'ajustement et en utilisant la courbe ajustée pour prédire la valeur du point qui a été laissé. Répétez cette opération pour tous les points afin que chacun soit laissé de côté une fois. Ensuite, calculez l'erreur de validation de votre courbe finale (courbe équipée de tous les points) en tant que moyenne des erreurs de prédiction.

Cela ne vous indiquera que la sensibilité de votre modèle pour tout nouveau point de données. Par exemple, il ne vous dira pas à quel point votre modèle énergétique est inexact. Cependant, il s'agira d'une estimation d'erreur beaucoup plus réaliste, simple erreur d'ajustement.

Vous pouvez également tracer des erreurs de prédiction en fonction du volume si vous le souhaitez.

Jman
la source