%% Parameter Estimation Example
%% Kurt Bryan and SIMIODE
% A very simple example of fitting a function or model to data by using least squares.
% The Data: Here are some hypothetical data in the form of (t,y) pairs:

tvals = [1.1,1.9,2.3,4.1,5.5]
yvals = [1.24,0.83,0.71,0.29,0.15]

% A quick plot of (time, distance) pairs:

scatter(tvals,yvals)

% The Model and Sum of Squares: Let's fit a model f(t) = a*exp(b*t) to this data by adjusting a and b. First form a sum of squares

syms f(t,a,b)
f(t,a,b) = a*exp(b*t)
syms SS(a,b)
SS(a,b) = sum((f(tvals,a,b)-yvals).^2);

% The Optimal Choice for a and b: We want to minimize SS as a function of a and b. First, plot SS to get a sense of where the minimum might occur (or actually, we'll plot log(SS)). It's clear that
% 1<a and b<0.

fsurf(log(SS),[1 3 -1 0])

% Rotating the graph around shows a around 2 and b around -0.5 looks promising. So set d(SS)/da = 0 and d(SS)/db = 0 and use Matlab's vpasolve command to find a good solution

dSSdaeqn = diff(SS,a)==0;
dSSdbeqn = diff(SS,b)==0;
abbest = vpasolve([dSSdaeqn,dSSdbeqn],[a,b],[2; -0.5]) %Initial guess a = 2, b = -0.5

% The residual is

abest = abbest.a
bbest = abbest.b
SS(abest,bbest)

% Use these values in f(t) to plot and compare to the data

fplot(f(t,abest,bbest),[1.1 5.5],'-r')
hold on;
scatter(tvals,yvals);
hold off;

% Alternatively, we can minimize SS with respect to a and b by using Matlab's built-in "fminsearch" minimization algorithm.

SSf = matlabFunction(SS,'Vars',{[a,b]});  %Converts the symbolic function "SS" to a traditional matlab function "SSf" that takes vector argument [k P]
[ab,fval] = fminsearch(SSf,[2.0, -0.5]) %Initial guess a = 2.0, b = -0.5