Stratégie d'ajustement d'une fonction hautement non linéaire

12

Pour analyser les données d'une expérience de biophysique, j'essaie actuellement de faire un ajustement de courbe avec un modèle hautement non linéaire. La fonction modèle ressemble essentiellement à:

y=ax+bx1/2

Ici, la valeur de est particulièrement intéressante.b

Un tracé pour cette fonction:

Graphique de fonction

(Notez que la fonction de modèle est basée sur une description mathématique approfondie du système et semble fonctionner très bien --- c'est juste que les ajustements automatisés sont délicats).

Bien sûr, la fonction du modèle est problématique: les stratégies d'ajustement que j'ai essayées jusqu'à présent échouent en raison de l'asymptote nette à , en particulier avec des données bruyantes.x=0

Ma compréhension du problème ici est que l'ajustement des moindres carrés simples (j'ai joué avec une régression linéaire et non linéaire dans MATLAB; principalement Levenberg-Marquardt) est très sensible à l'asymptote verticale, car les petites erreurs dans x sont énormément amplifiées .

Quelqu'un pourrait-il m'indiquer une stratégie appropriée qui pourrait contourner cela?

J'ai quelques connaissances de base en statistiques, mais c'est encore assez limité. Je serais impatient d'apprendre, si seulement je savais où commencer à chercher :)

Merci beaucoup pour vos conseils!

Modifier Je vous demande pardon d'avoir oublié de mentionner les erreurs. Le seul bruit significatif est en , et c'est additif.x

Edit 2 Quelques informations supplémentaires sur l'arrière-plan de cette question. Le graphique ci-dessus modélise le comportement d'étirement d'un polymère. Comme @whuber l'a souligné dans les commentaires, vous avez besoin de pour obtenir un graphique comme ci-dessus.b200a

Quant à la façon dont les gens ont ajusté cette courbe jusqu'à ce point: il semble que les gens coupent généralement l'asymptote verticale jusqu'à ce qu'ils trouvent un bon ajustement. Le choix de la coupure reste cependant arbitraire, ce qui rend la procédure d'ajustement peu fiable et non reproductible.

Modifier le graphique fixe 3 & 4 .

onnodb
la source
3
Les erreurs viennent-elles en ou en ou les deux? Sous quelle forme vous attendez-vous à ce que le bruit entre (multiplicatif, additif, etc.)? yxy
probabilités
2
@onnodb: Ma préoccupation est la suivante: cela ne remet-il pas fondamentalement en cause la robustesse de votre modèle lui-même? Quelle que soit la stratégie d'ajustement que vous utilisez, ne restera- t-il pas très sensible? Pouvez-vous avoir une confiance élevée dans une telle estimation de ? bbb
curious_cat
1
Malheureusement, cela ne fonctionnera toujours pas. Il n'y a tout simplement pas de combinaison possible de et qui reproduira même qualitativement le graphique que vous avez dessiné. ( Il est évident est négatif. doit être inférieur au moins la pente dans le graphique, mais positif, ce qui le met dans un intervalle étroit. Mais quand est dans cet intervalle, il est tout simplement pas assez grand pour surmonter l'énorme pic négatif à l'origine introduite par le terme .) Qu'avez-vous dessiné? Les données? Une autre fonction? b b a a b x 1 / deuxabbaabx1/2
whuber
1
Merci, mais c'est toujours faux. En étendant la tangente à ce graphique vers l'arrière à partir de tout point où , vous intercepterez l'axe y à . Étant donné que la pointe descendante à montre que est négatif, cette ordonnée à l'origine doit également être négative. Mais dans votre figure, il est parfaitement clair que la plupart de ces interceptions sont positives, allant jusqu'à . Ainsi , il est mathématiquement impossible qu'une équation comme peut décrire votre courbe , pas même environ. Au minimum, vous devez ajuster quelque chose comme .x > 0 ( 0 , 3 b / ( 2 x 1 / 2 ) ) 0 b 15,5 y = a x + b x 1 / 2 y = a x + b x 1 / 2 + c(x,ax+bx1/2)x>0(0,3b/(2x1/2))0b15.5y=ax+bx1/2y=ax+bx1/2+c
whuber
1
Avant de travailler sur cela, je voulais m'assurer de l'énoncé de la question: c'est pourquoi il est important d'avoir la fonction correcte. Je n'ai pas le temps de donner une réponse complète maintenant, mais je voudrais faire remarquer que "d'autres personnes" peuvent se tromper - mais cela dépend encore plus de détails, hélas. Si votre erreur est vraiment additive, il me semble qu'elle doit encore être fortement hétéroscédastique, sinon sa variance à de petites valeurs de serait vraiment minuscule. Que pouvez-vous nous dire - quantitativement - à propos de cette erreur? xxx
whuber

Réponses:

10

Les méthodes que nous utiliserions pour ajuster cela manuellement (c'est-à-dire l'analyse exploratoire des données) peuvent fonctionner remarquablement bien avec de telles données.

Je souhaite re- paramétrer légèrement le modèle afin de rendre ses paramètres positifs:

y=axb/x.

Pour un donné , supposons qu'il existe un vrai x unique satisfaisant cette équation; appelons cela f ( y ; a , b ) ou, pour être bref, f ( y ) lorsque ( a , b ) sont compris.yxf(y;a,b)f(y)(a,b)

Nous observons une collection de paires ordonnées où les x i s'écartent de f ( y i ; a , b ) par des variations aléatoires indépendantes avec des moyennes nulles. Dans cette discussion, je suppose qu'ils ont tous une variance commune, mais une extension de ces résultats (en utilisant les moindres carrés pondérés) est possible, évidente et facile à mettre en œuvre. Voici un exemple simulé d'une telle collection de 100 valeurs, avec a = 0,0001 , b = 0,1 , et une variance commune de σ(xi,yi)xif(yi;a,b)100a=0.0001b=0.1 .σ2=4

Tracé de données

Il s'agit d'un exemple (délibérément) difficile, comme peuvent l'apprécier les valeurs non physiques (négatives) et leur extraordinaire dispersion (qui est généralement de ± 2 unités horizontales , mais peut aller jusqu'à 5 ou 6 sur l' axe x ). Si nous pouvons obtenir un ajustement raisonnable à ces données qui se rapproche de l'estimation des a , b et σ 2 utilisés, nous aurons bien réussi.x±2 56xabσ2

Un raccord exploratoire est itératif. Chaque étage est constitué de deux étapes: estimer (basé sur les données et les estimations précédentes un et b de a et b , à partir de laquelle les valeurs précédentes prédites x Au i peut être obtenue pour le x i ) et ensuite estimer b . Étant donné que les erreurs sont en x , les ajustements estiment le x i à partir de ( y i ) , plutôt que l'inverse. Ordonner d'abord les erreurs dans x , quand xaa^b^abx^ixibxi(yi)xx est suffisamment grand,

xi1a(yi+b^x^i).

Par conséquent, nous pouvons mettre à jour un en ajustant ce modèle avec des moindres carrés (avis , il n'a qu'un seul paramètre - une pente, une --et pas ordonnée à l' origine) et en prenant l'inverse du coefficient que l'estimation mise à jour d' un .a^aa

Ensuite, lorsque est suffisamment petit, le terme quadratique inverse domine et nous trouvons (là encore au premier ordre dans les erreurs) quex

xib212a^b^x^3/2yi2.

Encore une fois des moindres carrés (avec juste un terme de pente ) on obtient une estimation mise à jour b par la racine carrée de la pente aménagée.bb^

Pour voir pourquoi cela fonctionne, une approximation exploratoire grossière de cet ajustement peut être obtenue en traçant contre 1 / y 2 i pour le plus petit x i . Mieux encore, parce que les x i sont mesurés avec erreur et les y i changent de façon monotone avec les x i , nous devons nous concentrer sur les données avec les plus grandes valeurs de 1 / y 2 i . Voici un exemple de notre jeu de données simulé montrant la plus grande moitié du y ixi1/yi2xixiyixi1/yi2yi en rouge, la plus petite moitié en bleu, et une ligne passant par l'origine correspond aux points rouges.

Figure

Les points s'alignent approximativement, bien qu'il y ait un peu de courbure aux petites valeurs de et y . (Remarquez le choix des axes: parce que x est la mesure, il est classique de la tracer sur l' axe vertical .) En focalisant l'ajustement sur les points rouges, où la courbure doit être minimale, nous devons obtenir une estimation raisonnable de b . La valeur de 0,096 indiquée dans le titre est la racine carrée de la pente de cette ligne: elle n'est que de 4 % inférieure à la valeur réelle!xyxb0.0964

À ce stade, les valeurs prévues peuvent être mises à jour via

x^i=f(yi;a^,b^).

Itérer jusqu'à ce que les estimations se stabilisent (ce qui n'est pas garanti) ou qu'elles parcourent de petites plages de valeurs (qui ne peuvent toujours pas être garanties).

axba^=0.0001960.0001b^=0.10730.1). Ce graphique montre une fois de plus les données sur lesquelles sont superposées (a) la vraie courbe en gris (en pointillés) et (b) la courbe estimée en rouge (solide):

S'adapte

3.734

Il y a quelques problèmes avec cette approche:

  • Les estimations sont biaisées. Le biais devient apparent lorsque l'ensemble de données est petit et que relativement peu de valeurs sont proches de l'axe des x. L'ajustement est systématiquement un peu bas.

  • yiyi

  • ab


Code

Ce qui suit est écrit en Mathematica .

estimate[{a_, b_, xHat_}, {x_, y_}] := 
  Module[{n = Length[x], k0, k1, yLarge, xLarge, xHatLarge, ySmall, 
    xSmall, xHatSmall, a1, b1, xHat1, u, fr},
   fr[y_, {a_, b_}] := Root[-b^2 + y^2 #1 - 2 a y #1^2 + a^2 #1^3 &, 1];
   k0 = Floor[1 n/3]; k1 = Ceiling[2 n/3];(* The tuning constants *)
   yLarge = y[[k1 + 1 ;;]]; xLarge = x[[k1 + 1 ;;]]; xHatLarge = xHat[[k1 + 1 ;;]];
   ySmall = y[[;; k0]]; xSmall = x[[;; k0]]; xHatSmall = xHat[[;; k0]];
   a1 = 1/
     Last[LinearModelFit[{yLarge + b/Sqrt[xHatLarge], 
          xLarge}\[Transpose], u, u]["BestFitParameters"]];
   b1 = Sqrt[
     Last[LinearModelFit[{(1 - 2 a1 b  xHatSmall^(3/2)) / ySmall^2, 
          xSmall}\[Transpose], u, u]["BestFitParameters"]]];
   xHat1 = fr[#, {a1, b1}] & /@ y;
   {a1, b1, xHat1}
   ];

xydata = {x,y}a=b=0

{a, b, xHat} = NestWhile[estimate[##, data] &, {0, 0, data[[1]]}, 
                Norm[Most[#1] - Most[#2]] >= 0.001 &,  2, 100]
whuber
la source
3
Ceci est une réponse étonnante; Je suis bien obligé! J'ai joué avec ça, et les résultats semblent très prometteurs. J'aurai besoin d'un peu plus de temps pour bien comprendre le raisonnement, cependant :) Aussi: pourrais-je vous contacter via votre site Web pour une question supplémentaire (privée), concernant les remerciements?
onnodb
3

Voir les questions importantes @probabilityislogic postées

y=yxyx=x3/21/x

b

x

-

Modifiez pour prendre en compte les informations supplémentaires:

y=b+ax

Nous avons maintenant que les erreurs sont en x et additives. Nous ne savons toujours pas si la variance est constante à cette échelle.

x=y/ab/a=my+c

xo=x+ηx

oxo

xo=c+my+ϵϵ=ζxy

Je ne suis pas sûr que cela améliore les choses! Je crois qu'il existe des méthodes pour ce genre de chose, mais ce n'est pas du tout mon domaine.

J'ai mentionné dans les commentaires que vous aimeriez peut-être examiner la régression inverse, mais la forme particulière de votre fonction peut vous empêcher d'aller trop loin.

Vous pourriez même être coincé avec des méthodes assez robustes à erreurs dans x sous cette forme linéaire.

-

y

x

Glen_b -Reinstate Monica
la source
x
2
" même si les erreurs sont en x " - yikes, c'est assez important. Vous voudrez peut-être vérifier la régression inverse.
Glen_b -Reinstate Monica
3
x=13(2ya+21/3y2(27a4b22a3y3+3327a8b44a7b2y3)1/3+(27a4b22a3y3+3327a8b44a7b2y3)1/321/3a2)
xoxox+ζx=(thatmonster)+ϵϵ=ζ
x(y)yb
0

Après quelques semaines d'expérimentation, une technique différente semble fonctionner le mieux dans ce cas particulier: l' ajustement des moindres carrés totaux . C'est une variante de l'ajustement habituel (non linéaire) des moindres carrés, mais au lieu de mesurer les erreurs d'ajustement le long d'un seul des axes (ce qui pose des problèmes dans des cas hautement non linéaires comme celui-ci), il prend en compte les deux axes.

Il existe une pléthore d'articles, de tutoriels et de livres disponibles sur le sujet, bien que le cas non linéaire soit plus difficile à cerner. Il y a même du code MATLAB disponible.

onnodb
la source
yy
@whuber Merci d'avoir exprimé vos préoccupations! En ce moment, je travaille toujours sur l'exécution de simulations pour sonder la fiabilité de l'adaptation TLS à ce problème. Ce que j'ai vu jusqu'à présent, cependant, c'est que la prise en compte par TLS des deux variables aide grandement à surmonter la non-linéarité élevée du modèle. Les ajustements de données simulées sont fiables et convergent très bien. Plus de travail doit être fait cependant, et je devrai certainement empiler votre méthode jusqu'à celle-ci, une fois que nous aurons plus de données réelles disponibles --- et examiner en détail vos préoccupations.
onnodb
OK - n'oubliez pas que j'ai des inquiétudes comparables sur la méthode que j'ai proposée!
whuber