Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Ajuste de datos no lineal

Este ejemplo muestra cómo ajustar una función no lineal a los datos utilizando varios algoritmos de Optimization Toolbox™.

Problema de configuración

Tenga en cuenta los siguientes datos:

Data = ...   [0.0000    5.8955    0.1000    3.5639    0.2000    2.5173    0.3000    1.9790    0.4000    1.8990    0.5000    1.3938    0.6000    1.1359    0.7000    1.0096    0.8000    1.0343    0.9000    0.8435    1.0000    0.6856    1.1000    0.6100    1.2000    0.5392    1.3000    0.3946    1.4000    0.3903    1.5000    0.5474    1.6000    0.3459    1.7000    0.1370    1.8000    0.2211    1.9000    0.1704    2.0000    0.2636];

Vamos a trazar estos puntos de datos.

t = Data(:,1); y = Data(:,2); % axis([0 2 -0.5 6]) % hold on plot(t,y,'ro') title('Data points')

% hold off

Nos gustaría encajar la función

y = c(1)*exp(-lam(1)*t) + c(2)*exp(-lam(2)*t)

a los datos.

Enfoque de solución mediantelsqcurvefit

La función resuelve este tipo de problema fácilmente.lsqcurvefit

Para empezar, defina los parámetros en términos de una variable x:

x(1) = c(1)

x(2) = lam(1)

x(3) = c(2)

x(4) = lam(2)

A continuación, defina la curva como una función de los parámetros x y los datos t:

F = @(x,xdata)x(1)*exp(-x(2)*xdata) + x(3)*exp(-x(4)*xdata);

Hemos establecido arbitrariamente nuestro punto inicial x0 de la siguiente manera: c (1) = 1, Lam (1) = 1, c (2) = 1, Lam (2) = 0:

x0 = [1 1 1 0];

Ejecutamos el solucionador y trazamos el ajuste resultante.

[x,resnorm,~,exitflag,output] = lsqcurvefit(F,x0,t,y)
Local minimum possible.  lsqcurvefit stopped because the final change in the sum of squares relative to  its initial value is less than the value of the function tolerance. 
x = 1×4

    3.0068   10.5869    2.8891    1.4003

resnorm = 0.1477 
exitflag = 3 
output = struct with fields:
    firstorderopt: 7.8881e-06
       iterations: 6
        funcCount: 35
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 0.0096
          message: '...'

 hold on plot(t,F(x,t)) hold off

Enfoque de solución mediantefminunc

Para resolver el problema utilizando, se establece la función objetivo como la suma de los cuadrados de los residuos.fminunc

Fsumsquares = @(x)sum((F(x,t) - y).^2); opts = optimoptions('fminunc','Algorithm','quasi-newton'); [xunc,ressquared,eflag,outputu] = ...     fminunc(Fsumsquares,x0,opts)
Local minimum found.  Optimization completed because the size of the gradient is less than the value of the optimality tolerance. 
xunc = 1×4

    2.8890    1.4003    3.0069   10.5862

ressquared = 0.1477 
eflag = 1 
outputu = struct with fields:
       iterations: 30
        funcCount: 185
         stepsize: 0.0017
     lssteplength: 1
    firstorderopt: 2.9476e-05
        algorithm: 'quasi-newton'
          message: '...'

Observe que encontró la misma solución que, pero tomó muchas más evaluaciones de funciones para hacerlo.fminunclsqcurvefit Los parámetros para están en el orden opuesto a los de; la mayor Lam es Lam (2), no Lam (1).fminunclsqcurvefit Esto no es sorprendente, el orden de las variables es arbitrario.

fprintf(['There were %d iterations using fminunc,' ...     ' and %d using lsqcurvefit.\n'], ...     outputu.iterations,output.iterations)
There were 30 iterations using fminunc, and 6 using lsqcurvefit. 
fprintf(['There were %d function evaluations using fminunc,' ...     ' and %d using lsqcurvefit.'], ...     outputu.funcCount,output.funcCount)
There were 185 function evaluations using fminunc, and 35 using lsqcurvefit. 

División de los problemas lineales y no lineales

Observe que el problema de ajuste es lineal en los parámetros c (1) y c (2). Esto significa que para cualquier valor de Lam (1) y Lam (2), podemos utilizar el operador de barra diagonal inversa para encontrar los valores de c (1) y c (2) que resuelven el problema de mínimos cuadrados.

Ahora rehacemos el problema como un problema bidimensional, buscando los mejores valores de Lam (1) y Lam (2). Los valores de c (1) y c (2) se calculan en cada paso utilizando el operador de barra diagonal inversa como se describió anteriormente.

type fitvector
function yEst = fitvector(lam,xdata,ydata) %FITVECTOR  Used by DATDEMO to return value of fitting function. %   yEst = FITVECTOR(lam,xdata) returns the value of the fitting function, y %   (defined below), at the data points xdata with parameters set to lam. %   yEst is returned as a N-by-1 column vector, where N is the number of %   data points. % %   FITVECTOR assumes the fitting function, y, takes the form % %     y =  c(1)*exp(-lam(1)*t) + ... + c(n)*exp(-lam(n)*t) % %   with n linear parameters c, and n nonlinear parameters lam. % %   To solve for the linear parameters c, we build a matrix A %   where the j-th column of A is exp(-lam(j)*xdata) (xdata is a vector). %   Then we solve A*c = ydata for the linear least-squares solution c, %   where ydata is the observed values of y.  A = zeros(length(xdata),length(lam));  % build A matrix for j = 1:length(lam)    A(:,j) = exp(-lam(j)*xdata); end c = A\ydata; % solve A*c = y for linear parameters c yEst = A*c; % return the estimated response based on c 

Resolver el problema utilizando, a partir de un punto inicial bidimensional Lam (1), Lam (2):lsqcurvefit

x02 = [1 0]; F2 = @(x,t) fitvector(x,t,y);  [x2,resnorm2,~,exitflag2,output2] = lsqcurvefit(F2,x02,t,y)
Local minimum possible.  lsqcurvefit stopped because the final change in the sum of squares relative to  its initial value is less than the value of the function tolerance. 
x2 = 1×2

   10.5861    1.4003

resnorm2 = 0.1477 
exitflag2 = 3 
output2 = struct with fields:
    firstorderopt: 4.4107e-06
       iterations: 10
        funcCount: 33
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 0.0080
          message: '...'

La eficiencia de la solución bidimensional es similar a la de la solución de cuatro dimensiones:

fprintf(['There were %d function evaluations using the 2-d ' ...     'formulation, and %d using the 4-d formulation.'], ...     output2.funcCount,output.funcCount)
There were 33 function evaluations using the 2-d formulation, and 35 using the 4-d formulation. 

El problema de división es más robusto a la suposición inicial

La elección de un punto de partida erróneo para el problema de cuatro parámetros original conduce a una solución local que no es global. La elección de un punto de partida con los mismos valores de Lam (1) y Lam (2) defectuosos para el problema de dos parámetros divididos conduce a la solución global. Para mostrar esto, reejecutamos el problema original con un punto de inicio que conduce a una solución local relativamente mala y comparamos el ajuste resultante con la solución global.

x0bad = [5 1 1 0]; [xbad,resnormbad,~,exitflagbad,outputbad] = ...     lsqcurvefit(F,x0bad,t,y)
Local minimum possible.  lsqcurvefit stopped because the final change in the sum of squares relative to  its initial value is less than the value of the function tolerance. 
xbad = 1×4

  -20.2519    2.4796   25.3757    2.4795

resnormbad = 2.2173 
exitflagbad = 3 
outputbad = struct with fields:
    firstorderopt: 0.0036
       iterations: 31
        funcCount: 160
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 0.0012
          message: '...'

 hold on plot(t,F(xbad,t),'g') legend('Data','Global fit','Bad local fit','Location','NE') hold off

 fprintf(['The residual norm at the good ending point is %f,' ...    ' and the residual norm at the bad ending point is %f.'], ...    resnorm,resnormbad)
The residual norm at the good ending point is 0.147723, and the residual norm at the bad ending point is 2.217300. 

Temas relacionados