Mínims quadrats no lineals

Els mínims quadrats no lineals són la forma d'anàlisi de mínims quadrats que s'utilitza per ajustar un conjunt de m observacions amb un model que és no lineal en n paràmetres desconeguts (m ≥ n). S'utilitza en algunes formes de regressió no lineal. La base del mètode és aproximar el model mitjançant un model lineal i refinar els paràmetres mitjançant iteracions successives. Hi ha moltes similituds amb els mínims quadrats lineals, però també algunes diferències significatives. En teoria econòmica, el mètode dels mínims quadrats no lineals s'aplica en (i) la regressió probit, (ii) la regressió llindar, (iii) la regressió suau, (iv) la regressió d'enllaç logístic, (v) els regressors transformats de Box-Cox (${\displaystyle m(x,\theta _{i})=\theta _{1}+\theta _{2}x^{(\theta _{3})}}$).^[1]

Considereu un conjunt de ${\displaystyle m}$ punts de dades, ${\displaystyle (x_{1},y_{1}),(x_{2},y_{2}),\dots ,(x_{m},y_{m}),}$ i una corba (funció model) ${\displaystyle {\hat {y}}=f(x,{\boldsymbol {\beta }}),}$ que a més de la variable ${\displaystyle x}$ també depèn de ${\displaystyle n}$ paràmetres, ${\displaystyle {\boldsymbol {\beta }}=(\beta _{1},\beta _{2},\dots ,\beta _{n}),}$ amb ${\displaystyle m\geq n.}$ Es vol trobar el vector ${\displaystyle {\boldsymbol {\beta }}}$ de paràmetres tals que la corba s'ajusti millor a les dades donades en el sentit dels mínims quadrats, és a dir, la suma dels quadrats^[2]$${\displaystyle S=\sum _{i=1}^{m}r_{i}^{2}}$$ es minimitza, on els residuals (errors de predicció dins la mostra) r_i es donen per $${\displaystyle r_{i}=y_{i}-f(x_{i},{\boldsymbol {\beta }})}$$ per ${\displaystyle i=1,2,\dots ,m.}$

El valor mínim de S es produeix quan el gradient és zero. Com que el model conté n paràmetres, hi ha n equacions de gradient: $${\displaystyle {\frac {\partial S}{\partial \beta _{j}}}=2\sum _{i}r_{i}{\frac {\partial r_{i}}{\partial \beta _{j}}}=0\quad (j=1,\ldots ,n).}$$En un sistema no lineal, les derivades

${\textstyle {\frac {\partial r_{i}}{\partial \beta _{j}}}}$ són funcions tant de la variable independent com dels paràmetres, de manera que en general aquestes equacions de gradient no tenen una solució tancada. En canvi, s'han de triar valors inicials per als paràmetres. A continuació, els paràmetres es refinen iterativament, és a dir, els valors s'obtenen per aproximació successiva, $${\displaystyle \beta _{j}\approx \beta _{j}^{k+1}=\beta _{j}^{k}+\Delta \beta _{j}.}$$

Aquí, k és un nombre d'iteració i el vector d'increments, ${\displaystyle \Delta {\boldsymbol {\beta }}}$ es coneix com el vector de desplaçament. A cada iteració, el model es linealitza per aproximació a una expansió polinòmica de Taylor de primer ordre al voltant de ${\displaystyle {\boldsymbol {\beta }}^{k}}$$${\displaystyle f(x_{i},{\boldsymbol {\beta }})\approx f(x_{i},{\boldsymbol {\beta }}^{k})+\sum _{j}{\frac {\partial f(x_{i},{\boldsymbol {\beta }}^{k})}{\partial \beta _{j}}}\left(\beta _{j}-\beta _{j}^{k}\right)=f(x_{i},{\boldsymbol {\beta }}^{k})+\sum _{j}J_{ij}\,\Delta \beta _{j}.}$$ La matriu jacobiana, J, és una funció de constants, la variable independent i els paràmetres, de manera que canvia d'una iteració a la següent. Així, en termes del model linealitzat,

$${\displaystyle {\frac {\partial r_{i}}{\partial \beta _{j}}}=-J_{ij}}$$ i els residus vénen donats per $${\displaystyle \Delta y_{i}=y_{i}-f(x_{i},{\boldsymbol {\beta }}^{k}),}$$$${\displaystyle r_{i}=y_{i}-f(x_{i},{\boldsymbol {\beta }})=\left(y_{i}-f(x_{i},{\boldsymbol {\beta }}^{k})\right)+\left(f(x_{i},{\boldsymbol {\beta }}^{k})-f(x_{i},{\boldsymbol {\beta }})\right)\approx \Delta y_{i}-\sum _{s=1}^{n}J_{is}\Delta \beta _{s}.}$$

Si substituïm aquestes expressions a les equacions de gradient, esdevenen $${\displaystyle -2\sum _{i=1}^{m}J_{ij}\left(\Delta y_{i}-\sum _{s=1}^{n}J_{is}\ \Delta \beta _{s}\right)=0,}$$ que, en reorganitzar-se, es converteixen n equacions lineals simultànies, les equacions normals

$${\displaystyle \sum _{i=1}^{m}\sum _{s=1}^{n}J_{ij}J_{is}\ \Delta \beta _{s}=\sum _{i=1}^{m}J_{ij}\ \Delta y_{i}\qquad (j=1,\dots ,n).}$$Les equacions normals s'escriuen en notació matricial com a

$${\displaystyle \left(\mathbf {J} ^{\mathsf {T}}\mathbf {J} \right)\Delta {\boldsymbol {\beta }}=\mathbf {J} ^{\mathsf {T}}\ \Delta \mathbf {y} .}$$Aquestes equacions formen la base de l'algoritme de Gauss-Newton per a un problema de mínims quadrats no lineal.

Tingueu en compte la convenció de signes en la definició de la matriu jacobiana en termes de les derivades. Fórmules lineals en ${\displaystyle J}$ pot aparèixer amb un factor de ${\displaystyle -1}$ en altres articles o en la literatura.

En els mínims quadrats lineals, la funció objectiu, S, és una funció quadràtica dels paràmetres. $${\displaystyle S=\sum _{i}W_{ii}\left(y_{i}-\sum _{j}X_{ij}\beta _{j}\right)^{2}}$$ Quan només hi ha un paràmetre, el gràfic de S respecte a aquest paràmetre serà una paràbola. Amb dos o més paràmetres, els contorns de S respecte a qualsevol parell de paràmetres seran el·lipses concèntriques (suposant que la matriu de les equacions normals ${\displaystyle \mathbf {X} ^{\mathsf {T}}\mathbf {WX} }$ és definit positiu ). Els valors mínims dels paràmetres es troben al centre de les el·lipses. La geometria de la funció objectiu general es pot descriure com a el·líptica paraboloide. En NLLSQ la funció objectiu és quadràtica respecte als paràmetres només en una regió propera al seu valor mínim, on la sèrie de Taylor truncada és una bona aproximació al model. $${\displaystyle S\approx \sum _{i}W_{ii}\left(y_{i}-\sum _{j}J_{ij}\beta _{j}\right)^{2}}$$ Com més difereixen els valors dels paràmetres dels seus valors òptims, més es desvien els contorns de la forma el·líptica. Una conseqüència d'això és que les estimacions inicials dels paràmetres han de ser el més properes possible als seus valors òptims (desconeguts!). També explica com es pot produir la divergència, ja que l'algoritme de Gauss-Newton només és convergent quan la funció objectiu és aproximadament quadràtica en els paràmetres.

Les equacions normals $${\displaystyle \left(\mathbf {J} ^{\mathsf {T}}\mathbf {WJ} \right)\Delta {\boldsymbol {\beta }}=\left(\mathbf {J} ^{\mathsf {T}}\mathbf {W} \right)\Delta \mathbf {y} }$$ es pot resoldre per a ${\displaystyle \Delta {\boldsymbol {\beta }}}$ per descomposició de Cholesky, tal com es descriu en mínims quadrats lineals. Els paràmetres s'actualitzen iterativament. $${\displaystyle {\boldsymbol {\beta }}^{k+1}={\boldsymbol {\beta }}^{k}+\Delta {\boldsymbol {\beta }}}$$ on k és un nombre d'iteració. Tot i que aquest mètode pot ser adequat per a models simples, fallarà si es produeix divergència. Per tant, la protecció contra la divergència és essencial.

El mínim en la suma de quadrats es pot trobar mitjançant un mètode que no implica la formació de les equacions normals. Els residuals amb el model linealitzat es poden escriure com $${\displaystyle \mathbf {r} =\Delta \mathbf {y} -\mathbf {J} \,\Delta {\boldsymbol {\beta }}.}$$ El Jacobià se sotmet a una descomposició ortogonal; la descomposició QR servirà per il·lustrar el procés. $${\displaystyle \mathbf {J} =\mathbf {QR} }$$ on Q és ortogonal ${\displaystyle m\times m}$ matriu i R és una ${\displaystyle m\times n}$ matriu que es divideix en una ${\displaystyle n\times n}$ bloc, ${\displaystyle \mathbf {R} _{n}}$, i un ${\displaystyle (m-n)\times n}$ bloc zero. ${\displaystyle \mathbf {R} _{n}}$ és triangular superior.

Una variant del mètode de descomposició ortogonal implica la descomposició de valors singulars, en què R es diagonalitza mitjançant transformacions ortogonals addicionals.

Hi ha molts exemples a la literatura científica on s'han utilitzat diferents mètodes per a problemes d'ajust de dades no lineals.

• Inclusió de segones derivades a l'expansió en sèrie de Taylor de la funció model. Aquest és el mètode de Newton en optimització. $${\displaystyle f(x_{i},{\boldsymbol {\beta }})=f^{k}(x_{i},{\boldsymbol {\beta }})+\sum _{j}J_{ij}\,\Delta \beta _{j}+{\frac {1}{2}}\sum _{j}\sum _{k}\Delta \beta _{j}\,\Delta \beta _{k}\,H_{jk_{(i)}},\ H_{jk_{(i)}}={\frac {\partial ^{2}f(x_{i},{\boldsymbol {\beta }})}{\partial \beta _{j}\,\partial \beta _{k}}}.}$$ La matriu H es coneix com a matriu hessiana. Tot i que aquest model té millors propietats de convergència prop del mínim, és molt pitjor quan els paràmetres són lluny dels seus valors òptims. El càlcul de la matriu hessiana augmenta la complexitat de l'algoritme. Aquest mètode no és d'ús general.
• Mètode de Davidon-Fletcher-Powell. Aquest mètode, una forma del mètode pseudo-Newton, és similar a l'anterior però calcula el hessià per aproximació successiva, per evitar haver d'utilitzar expressions analítiques per a les segones derivades.
• Descens més pronunciat. Tot i que es garanteix una reducció de la suma de quadrats quan el vector de desplaçament apunta en la direcció del descens més pronunciat, aquest mètode sovint té un rendiment deficient. Quan els valors dels paràmetres estan lluny de ser òptims, la direcció del vector de descens més pronunciat, que és normal (perpendicular) als contorns de la funció objectiu, és molt diferent de la direcció del vector de Gauss-Newton. Això fa que la divergència sigui molt més probable, sobretot perquè el mínim al llarg de la direcció del descens més pronunciat pot correspondre a una petita fracció de la longitud del vector de descens més pronunciat. Quan els contorns de la funció objectiu són molt excèntrics, a causa d'una alta correlació entre els paràmetres, les iteracions de descens més pronunciat, amb tall per desplaçament, segueixen una trajectòria lenta i en ziga-zaga cap al mínim.
• Cerca de gradient conjugat. Aquest és un mètode millorat basat en el descens més pronunciat amb bones propietats de convergència teòrica, tot i que pot fallar en ordinadors digitals de precisió finita fins i tot quan s'utilitza en problemes quadràtics.

Els mètodes de cerca directa depenen d'avaluacions de la funció objectiu en una varietat de valors de paràmetres i no utilitzen derivades en absolut. Ofereixen alternatives a l'ús de derivades numèriques en el mètode de Gauss-Newton i els mètodes de gradient.

• Cerca de variables alternants. Cada paràmetre es varia al seu torn afegint-hi un increment fix o variable i mantenint el valor que provoca una reducció de la suma de quadrats. El mètode és senzill i eficaç quan els paràmetres no estan altament correlacionats. Té propietats de convergència molt deficients, però pot ser útil per trobar estimacions inicials de paràmetres.
• Cerca de Nelder-Mead (símplex). Un símplex en aquest context és un polítop de n + 1 vèrtexs en n dimensions; un triangle en un pla, un tetraedre en un espai tridimensional, etc. Cada vèrtex correspon a un valor de la funció objectiu per a un conjunt particular de paràmetres. La forma i la mida del simplex s'ajusten variant els paràmetres de manera que el valor de la funció objectiu al vèrtex més alt sempre disminueixi. Tot i que la suma de quadrats pot disminuir inicialment ràpidament, pot convergir a un punt no estacionari en problemes quasiconvexos, amb un exemple de MJD Powell.

Podeu trobar descripcions més detallades d'aquests i altres mètodes a Receptes numèriques, juntament amb codi informàtic en diversos idiomes.^[4]