Ajuste por mínimos cuadrados

La traducción de esta página aún no se ha actualizado a la versión más reciente. Haga clic aquí para ver la última versión en inglés.

Introducción

Curve Fitting Toolbox™ utiliza el método de mínimos cuadrados cuando ajusta datos. El ajuste requiere un modelo paramétrico que relacione los datos de respuesta con los datos predictores con uno o más coeficientes. El resultado del proceso de ajuste es una estimación de los coeficientes modelo.

Para obtener las estimaciones de coeficiente, el método de los mínimos cuadrados minimiza la suma de cuadrados de los valores residuales. El valor residual para el punto de datos i-ésimo r_i se define como la diferencia entre el valor de la respuesta observada y_i y el valor de la respuesta ajustada ŷ_i, y se identifica como error asociado con los datos.

$r_{i} = y_{i} - {\hat{y}}_{i}$

La suma de cuadrados de los valores residuales viene dada por

$S S E = \sum_{i = 1}^{n} r_{i}^{2} = \sum_{i = 1}^{n} {(y_{i} - {\hat{y}}_{i})}^{2}$

donde n es el número de puntos de datos incluidos en el ajuste y S es la suma de cuadrados de estimación de error. Los tipos compatibles para el ajuste por mínimos cuadrados incluyen:

Mínimos cuadrados lineales
Mínimos cuadrados lineales ponderados
Mínimos cuadrados robustos
Mínimos cuadrados no lineales

Distribuciones de error

Al ajustar datos que contienen variaciones aleatorias, hay dos supuestos importantes que se suelen hacer sobre el error:

El error existe solo en los datos de respuesta y no en los datos predictores.
Los errores son aleatorios y siguen una distribución normal (de Gauss) con media cero y varianza constante, σ².

El segundo supuesto se suele expresar como

$e r r o r \sim N (0, σ^{2})$

Se asume que los errores se distribuyen de forma normal porque la distribución normal suele proporcionar una aproximación adecuada a la distribución de muchas cantidades medidas. Aunque el método de ajuste por mínimos cuadrados no asume errores distribuidos con normalidad cuando calcula estimaciones de parámetro, el método funciona mejor para datos que no contienen un gran número de errores aleatorios con valores extremos. La distribución normal es una de las distribuciones de probabilidad en la que los errores aleatorios extremos no son habituales. Sin embargo, los resultados estadísticos como los límites de confianza y de predicción sí requieren errores distribuidos con normalidad para su validez.

Si la media de los errores es cero, entonces los errores son puramente aleatorios. Si la media no es cero, entonces puede que el modelo no sea la elección correcta para sus datos, o que los errores no sean puramente aleatorios y contengan errores sistemáticos.

Una varianza constante en los datos implica que la "difusión" de errores es constante. A veces se dice que los datos que tienen la misma varianza son de igual calidad.

El supuesto de que los errores aleatorios tienen una varianza constante no está implícito en la regresión de mínimos cuadrados ponderados. En cambio, se asume que las ponderaciones proporcionadas en el procedimiento de ajuste indican correctamente los distintos niveles de calidad presentes en los datos. Así, las ponderaciones se utilizan para ajustar a un nivel apropiado la influencia que tiene cada punto de datos en las estimaciones de los coeficientes ajustados.

Mínimos cuadrados lineales

Curve Fitting Toolbox utiliza el método de mínimos cuadrados para ajustar un modelo lineal a los datos. Un modelo lineal se define como una ecuación que es lineal en los coeficientes. Por ejemplo, los polinomiales son lineales pero los de Gauss no lo son. Para ilustrar el proceso de ajuste por mínimos cuadrados lineales, suponga que tiene n puntos de datos que pueden modelarse con un polinomio de primer grado.

$y = p_{1} x + p_{2}$

Para resolver esta ecuación para los coeficientes desconocidos p₁ y p₂, escriba S como un sistema de n ecuaciones lineales simultáneas en dos incógnitas. Si n es mayor que el número de incógnitas, entonces el sistema de ecuaciones está sobredeterminado.

$S S E = \sum_{i = 1}^{n} {(y_{i} - (p_{1} x_{i} + p_{2}))}^{2}$

Dado que el proceso de ajuste por mínimos cuadrados minimiza la suma de cuadrados de los valores residuales, los coeficientes se determinan diferenciando S con respecto a cada parámetro, y definiendo el resultado igual a cero.

$\begin{array}{l} \frac{\partial S S E}{\partial p_{1}} = - 2 \sum_{i = 1}^{n} x_{i} (y_{i} - (p_{1} x_{i} + p_{2})) = 0 \\ \frac{\partial S S E}{\partial p_{2}} = - 2 \sum_{i = 1}^{n} (y_{i} - (p_{1} x_{i} + p_{2})) = 0 \end{array}$

Las estimaciones de los parámetros verdaderos suelen representarse con b. Al sustituir b₁ y b₂ por p₁ y p₂, la anterior ecuación se convierte en

$\begin{array}{l} \sum_{i = 1}^{n} x_{i} (y_{i} - (b_{1} x_{i} + b_{2})) = 0 \\ \sum_{i = 1}^{n} (y_{i} - (b_{1} x_{i} + b_{2})) = 0 \end{array}$

en la que las sumas van de i = 1 para n. Las ecuaciones normales se definen como

$\begin{matrix} b_{1} \sum_{i = 1}^{n} x_{i}^{2} + b_{2} \sum_{i = 1}^{n} x_{i} = \sum_{i = 1}^{n} x_{i} y_{i} \\ b_{1} \sum_{i = 1}^{n} x_{i} + n b_{2} = \sum_{i = 1}^{n} y_{i} \end{matrix}$

Solución para b₁

$b_{1} = \frac{n \sum_{i = 1}^{n} x_{i} y_{i} - \sum_{i = 1}^{n} x_{i} \sum_{i = 1}^{n} y_{i}}{n \sum_{i = 1}^{n} x_{i}^{2} - {(\sum_{i = 1}^{n} x_{i})}^{2}}$

Solución para b₂ con el valor b₁

$b_{2} = \frac{1}{n} (\sum_{i = 1}^{n} y_{i} - b_{1} \sum_{i = 1}^{n} x_{i})$

Como puede observar, estimar los coeficientes p₁ y p₂ solo requiere unos pocos cálculos sencillos. Este ejemplo se puede extender a polinomios de mayor grado de forma sencilla aunque algo tediosa. Todo lo que se necesita es una ecuación normal adicional para cada término lineal que se añada al modelo.

En forma de matriz, los modelos lineales vienen dados por la fórmula

y = Xβ + ε

donde

y es un vector de respuestas n por 1.
β es un vector de coeficientes m por 1.
X es la matriz de diseño para el modelo n por m.
ε es un vector de errores n por 1.

Para el polinomio de primer grado, las ecuaciones n en dos incógnitas se expresan en términos de y, X y β como

$[\begin{matrix} y_{1} \\ y_{2} \\ y_{3} \\ . \\ . \\ . \\ y_{n} \end{matrix}] = [\begin{matrix} x_{1} 1 \\ x_{2} 1 \\ x_{3} 1 \\ . \\ . \\ . \\ x_{n} 1 \end{matrix}] \times [\begin{array}{l} p_{1} \\ p_{2} \end{array}]$

La solución de mínimos cuadrados al problema es un vector b que estima el vector desconocido de los coeficientes β. Las ecuaciones normales vienen dadas por

(X^TX)b = X^Ty

donde X^T es la traspuesta de la matriz de diseño X. Solución para b,

b = (X^TX)^–1 X^Ty

Utilice el operador de barra invertida de MATLAB^® (mldivide) para resolver un sistema de ecuaciones lineales simultáneas para coeficientes desconocidos. Dado que invertir X^TX puede conducir a errores de redondeo inaceptables, el operador de barra invertida utiliza la descomposición QR con pivoteo, que es un algoritmo muy estable numéricamente. Consulte Operaciones aritméticas para obtener más información sobre el operador de barra invertida y la descomposición QR.

Puede reincorporar b en la fórmula del modelo para obtener los valores de respuesta predichos, ŷ.

ŷ = Xb = Hy

H = X(X^TX)^–1 X^T

Un sombrero (acento circunflejo) sobre una letra denota una estimación de un parámetro o un pronóstico de un modelo. La matriz de proyección H se llama la matriz del sombrero porque coloca un sombrero sobre la y.

Los valores residuales vienen dados por

r = y – ŷ = (1–H)y

Mínimos cuadrados ponderados

Se suele asumir que los datos de respuesta son de igual calidad y, por lo tanto, tienen varianza constante. Si se infringe este supuesto, puede que esté indebidamente influenciado por datos de poca calidad. Para mejorar el ajuste, puede utilizar la regresión de mínimos cuadrados ponderados en la que se incluye un factor de escala adicional (la ponderación) en el proceso de ajuste. La regresión de mínimos cuadrados ponderados reduce la estimación de error

$S S E = \sum_{i = 1}^{n} w_{i} {(y_{i} - {\hat{y}}_{i})}^{2}$

donde w_i son las ponderaciones. Las ponderaciones determinan cuánto influye cada valor de respuesta en las estimaciones finales de los parámetros. Un punto de datos de gran calidad influye en el ajuste más que un punto de datos de baja calidad. Se recomienda ponderar los datos si se conocen las ponderaciones o si se justifica que siguen una forma concreta.

Las ponderaciones modifican la expresión para las estimaciones de parámetros b de la siguiente manera,

$b = {(X^{T} W X)}^{- 1} X^{T} W y$

donde W viene dada por los elementos diagonales de la matriz de ponderaciones w.

Ajustando los datos y representando los valores residuales, suele ser posible determinar si las varianzas no son constantes. En la gráfica que se muestra a continuación, los datos contienen datos replicados de distinta calidad y se asume que el ajuste es correcto. La gráfica de valores residuales revela los datos de poca calidad y tiene forma de "embudo", en el que los valores predictores pequeños ofrecen una mayor dispersión en los valores de respuesta que los valores predictores grandes.

Las ponderaciones que proporcione deberían transformar las varianzas de respuesta en un valor constante. Si conoce las varianzas de los errores de medición de sus datos, entonces las ponderaciones vienen dadas por

$w_{i} = 1 / σ_{i}^{2}$

Si solo dispone de las estimaciones de la variable de error para cada punto de datos, suele ser suficiente utilizar esas estimaciones en lugar de la verdadera varianza. Si no conoce las varianzas, basta con que especifique las ponderaciones en una escala relativa. Tenga en cuenta que un término de varianza general se estima incluso cuando se han especificado las ponderaciones. En este ejemplo, las ponderaciones definen la ponderación relativa de cada punto en el ajuste, pero no se toman para especificar la varianza exacta de cada punto.

Por ejemplo, si cada punto de datos es la media de varias mediciones independientes, podría tener sentido utilizar esos números de mediciones como ponderaciones.

Mínimos cuadrados robustos

Se suele asumir que los errores de respuesta siguen una distribución normal y que los valores extremos no son habituales. Aun así, los valores extremos llamados atípicos sí pueden aparecer.

La principal desventaja del ajuste por mínimos cuadrados es su sensibilidad a los valores atípicos. Los valores atípicos tienen una gran influencia en el ajuste porque cuadrar los valores residuales magnifica los efectos de estos puntos de datos extremos. Para minimizar la influencia de los valores atípicos, puede ajustar sus datos con la regresión de mínimos cuadrados robusta. La toolbox proporciona estos dos métodos de regresión robusta:

El mínimo de valores residuales absoluto (LAR). El método LAR identifica una curva que minimiza la diferencia absoluta de los valores residuales, en lugar de las diferencias cuadráticas. Por lo tanto, los valores extremos tienen una menor influencia en el ajuste.
Ponderaciones bicuadradas. Este método minimiza una suma ponderada de cuadrados, donde la ponderación dada a cada punto de datos depende de lo lejos que esté cada punto de la línea de ajuste. Los puntos cerca de la línea obtienen una ponderación completa. Los puntos más alejados de la línea obtienen una ponderación reducida. Los puntos que están más alejados de la línea de lo que cabría esperar por casualidad obtienen ponderación cero.
En la mayoría de los casos, el método de ponderación bicuadrada es preferible al LAR porque intenta buscar una curva que ajuste el grueso de los datos con el enfoque de mínimos cuadrados habitual y, al mismo tiempo, minimiza el efecto de los valores atípicos.

El ajuste robusto con ponderaciones bicuadradas utiliza un algoritmo de mínimos cuadrados reponderado de forma iterativa y sigue este procedimiento:

Ajusta el modelo por mínimos cuadrados ponderados.
Calcula los valores residuales ajustados y los estandariza. Los valores residuales ajustados vienen dados por
$r_{a d j} = \frac{r_{i}}{\sqrt{1 - h_{i}}}$
r_i son los valores residuales de mínimos cuadrados habituales y h_i son aprovechados que ajustan los valores residuales reduciendo la ponderación de los puntos de datos de alto aprovechamiento, que tienen un gran efecto en el ajuste por mínimos cuadrados. Los valores residuales ajustados estandarizados vienen dados por
$u = \frac{r_{a d j}}{K s}$
K es una constante de ajuste igual a 4,685 y s es la desviación robusta estándar dada por MAD/0,6745, donde MAD es la desviación absoluta mediana de los valores residuales.
Calcule las ponderaciones robustas como una función de u. Las ponderaciones bicuadradas vienen dadas por
$w_{i} = {\begin{matrix} {(1 - {(u_{i})}^{2})}^{2} & | u_{i} | < 1 \\ 0 & | u_{i} | \geq 1 \end{matrix}$
Tenga en cuenta que si proporciona su propio vector de ponderación regresiva, la ponderación será el producto de la ponderación robusta y la ponderación regresiva.
Si el ajuste converge, entonces habrá terminado. De lo contrario, puede realizar la siguiente iteración del procedimiento de ajuste regresando al primer paso.

La gráfica que se muestra a continuación compara un ajuste lineal regular con un ajuste robusto con ponderaciones bicuadradas. Tenga en cuenta que el ajuste robusto sigue el grueso de los datos y no está muy influenciado por los valores atípicos.

En lugar de minimizar los efectos de los valores atípicos con una regresión robusta, puede marcar puntos de datos que se excluyan del ajuste. Consulte Eliminar valores atípicos para obtener más información.

Mínimos cuadrados no lineales

El software de Curve Fitting Toolbox utiliza la formulación de mínimos cuadrados no lineales para ajustar un modelo no lineal a los datos. Un modelo no lineal se define como una ecuación que es no lineal en los coeficientes, o una combinación de lineales y no lineales en los coeficientes. Por ejemplo, las de Gauss, las relaciones de polinomios y las funciones de potencia son todas no lineales.

En forma de matriz, los modelos no lineales vienen dados por la fórmula

y = f(X,β) + ε

donde

y es un vector de respuestas n por 1.
f es una función de β y X.
β es un vector de coeficientes m por 1.
X es la matriz de diseño para el modelo n por m.
ε es un vector de errores n por 1.

Los modelos no lineales son más difíciles de ajustar que los modelos lineales porque los coeficientes no se pueden estimar con simples técnicas de matriz. En su lugar, se requiere un enfoque iterativo que siga estos pasos:

Comience con una estimación inicial para cada coeficiente. Para algunos modelos no lineales, se proporciona un enfoque heurístico que produce valores iniciales razonables. Para otros modelos, se proporcionan valores aleatorios del intervalo [0,1].
Produzca la curva de ajuste para el conjunto actual de coeficientes. El valor de respuesta ajustada ŷ viene dado por
ŷ = f(X,b)
e implica el cálculo de la jacobiana de f(X,b), que se define como la matriz de derivativas parciales tomadas con respecto a los coeficientes.
Ajuste los coeficientes y determine si el ajuste mejora. La dirección y la magnitud del ajuste dependen del algoritmo de ajuste. La toolbox proporciona estos algoritmos:
- Trust-region: es el algoritmo predeterminado y debe usarse si se especifican limitaciones de coeficiente. Puede resolver problemas no lineales con más eficiencia que el resto de algoritmos y representa una mejora sobre el popular algoritmo Levenberg-Marquardt.
- Levenberg-Marquardt: este algoritmo se ha usado durante años y ha demostrado que funcionaba la mayoría de las veces para una amplia gama de modelos no lineales y valores iniciales. Si el algoritmo de región de confianza no produce un ajuste razonable, y no dispone de las limitaciones de coeficiente, debería probar el algoritmo Levenberg-Marquardt.
Repita el proceso regresando al paso 2 hasta que el ajuste alcance los criterios de convergencia especificados.

Puede utilizar ponderaciones y ajuste robusto para modelos no lineales, y el proceso de ajuste se modificará en consecuencia.

Dada la naturaleza del proceso de aproximación, ningún algoritmo es infalible para todos los modelos no lineales, conjuntos de datos y puntos de partida. Por lo tanto, si no logra un ajuste razonable con los puntos de partida, algoritmo y criterio de convergencia predeterminados, debería probar con distintas opciones. Consulte Especificar las opciones de ajuste y los puntos de partida optimizados para obtener una descripción de cómo modificar las opciones predeterminadas. Dado que los modelos no lineales pueden ser especialmente sensibles a los puntos de partida, debería modificar esta opción en primer lugar.

Ajuste robusto

Abrir script en vivo

Este ejemplo muestra cómo comparar los efectos de excluir valores atípicos y ajuste robusto. El ejemplo muestra cómo excluir valores atípicos a una distancia arbitraria mayor que 1,5 desviaciones estándar del modelo. A continuación, los pasos comparan la acción de retirar los valores atípicos con la de especificar un ajuste robusto que ofrece ponderaciones menores para los valores atípicos.

Cree una señal sinusoidal de referencia:

xdata = (0:0.1:2*pi)'; 
y0 = sin(xdata);

Añada ruido a la señal con varianza no constante.

% Response-dependent Gaussian noise
gnoise = y0.*randn(size(y0));

% Salt-and-pepper noise
spnoise = zeros(size(y0)); 
p = randperm(length(y0));
sppoints = p(1:round(length(p)/5));
spnoise(sppoints) = 5*sign(y0(sppoints));

ydata = y0 + gnoise + spnoise;

Ajuste los datos ruidosos con un modelo sinusoidal de referencia y especifique tres argumentos de salida para obtener información de ajuste que incluya valores residuales.

f = fittype('a*sin(b*x)'); 
[fit1,gof,fitinfo] = fit(xdata,ydata,f,'StartPoint',[1 1]);

Analice la información en la estructura fitinfo.

fitinfo

fitinfo = struct with fields:
           numobs: 63
         numparam: 2
        residuals: [63x1 double]
         Jacobian: [63x2 double]
         exitflag: 3
    firstorderopt: 0.0883
       iterations: 5
        funcCount: 18
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 4.1539e-04
          message: 'Success, but fitting stopped because change in residuals less than tolerance (TolFun).'

Obtenga los valores residuales de la estructura fitinfo.

residuals = fitinfo.residuals;

Identifique los "valores atípicos" a una distancia arbitraria mayor que 1,5 desviaciones estándar del modelo de referencia y reajuste los datos con los valores atípicos excluidos.

I = abs( residuals) > 1.5 * std( residuals );
outliers = excludedata(xdata,ydata,'indices',I);

fit2 = fit(xdata,ydata,f,'StartPoint',[1 1],...
           'Exclude',outliers);

Compare el efecto de excluir los valores atípicos con el efecto de darles una ponderación bicuadrada menor en un ajuste robusto.

fit3 = fit(xdata,ydata,f,'StartPoint',[1 1],'Robust','on');

Represente los datos, los valores atípicos y los resultados de los ajustes. Especifique una leyenda informativa.

plot(fit1,'r-',xdata,ydata,'k.',outliers,'m*') 
hold on
plot(fit2,'c--')
plot(fit3,'b:')
xlim([0 2*pi])
legend( 'Data', 'Data excluded from second fit', 'Original fit',...
    'Fit with points excluded', 'Robust fit' )
hold off

Figure contains an axes object. The axes object contains 5 objects of type line. These objects represent Data, Data excluded from second fit, Original fit, Fit with points excluded, Robust fit.

Represente los valores residuales para los dos ajustes teniendo en cuenta los valores atípicos:

figure 
plot(fit2,xdata,ydata,'co','residuals') 
hold on
plot(fit3,xdata,ydata,'bx','residuals')
hold off

Figure contains an axes object. The axes object contains 4 objects of type line. These objects represent data, zero line.