Cómo ajustar un conjunto de datos a una función cualquiera

Tenemos un conjunto de puntos que queremos ajustar a una cierta función. Por ejemplo, tenemos los puntos:

x = [1 2 3 4 5]';
y = [1.24 0.59 0.31 0.19 0.09]';

Y queremos ajustarlos a la función y = k·exp(-a·x). Para ello, tenemos que estimar los parámetros "k" y "a".

---

Primero, creamos la función a ajustar en un archivo .m, por ejemplo "f01.m":

function z = f01(u,x)
z = u(1)*exp(-u(2)*x);

donde u(1) = k, y u(2) = a .

A continuación, creamos la función objetivo que vamos a optimizar, y que sume los errores al cuadro entre los valores medidos y los calculados. La siguiente función la creamos en el archivo "fobj.m":

function z = fobj(u,x,ym)

ycal = f01(u,x);
z = sum((ym-ycal).^2);

Por último, calculamos los parámetros "k" y "a" con la función "fminsearch":

[u,f] = fminsearch(@(u) fobj(u,x,y),[1 1])

donde [1 1] es el valor inicial de los parámetros, con los que el optimizador empezará a buscar la solución, la primera salida de la función "u" es un vector con los parámetros estimados, y la segunda salida de la función "f" es el valor de la función objetivo en el punto mínimo.

Matlab nos devuelve el siguiente resultado:

u =

    2.4430    0.6864

f =

    0.0022

---

En vez de usar "fminsearch" y tener que crear una función objetivo, podríamos haber usado directamente la funcion "lsqcurvefit". Pero la verdad, es que el método anterior me suele funcionar mejor.

[u,f] = lsqcurvefit(@(u,x) f01(u,x),[1,1],x,y)

donde [1 1] son los valores iniciales de u(1) y u(2), con los que la función empieza a tantear. Por cierto, cuando escribimos @(u,x), hay que poner primero el vector de parámetros a estimar, y segundo la variable independiente.

---

Por último, representamos tanto los puntos reales como la curva obtenida:

xx = [0:0.1:6];
yy = f01(u,xx);

figure(1)
hold on
box on
plot(xx,yy,'-','Color',[1 0.6 0.78],'LineWidth',2)
plot(x,y,'o','MarkerSize',5,'Color',[0.48 0.06 0.89],'MarkerFaceColor',[0.48 0.06 0.89])
legend({'real','estimado'},1,'EdgeColor',[1 1 1]);
axis([0 6 0 3])
set(gca,'XTick',[0:2:6])
set(gca,'YTick',[0:1:3])
set(gca,'Fontsize',10)

Nota: He usado 'EdgeColor',[1 1 1] para que la leyenda no esté encerrada en una caja. Y 'MarkerFaceColor',[0.48 0.06 0.89] para que los círculos estén rellenos.