How to use Octave to Fit an Arbitrary Function with fmins

From Class Wiki
Revision as of 17:38, 6 January 2011 by Frohro (talk | contribs)
Jump to navigation Jump to search

Suppose we want to fit experimental data to an arbitrary function. For our purposes let that function be . We can use the fmins function from the octave package optim. We are looking for a best fit and . The fmins function will try a whole lot of different values for these parameters until it decides to give up of it has found a local minimum.

The first step is to construct a function that computes the sum of the differences between the guess for the best fit function and the experimental data. This is the sum of the squared errors at each data point. The idea is that octave will use the fmins function to find the parameters that minimize this sum of squared errors. This can be but in a separate dot m file, or in the same script file you are using to do the fit.

function sum_square_errors = model(p,x,y)
% This is where we computer the sum of the square of the errors.
% The parameters are in the vector p, which for us is a two by one.
A=p(1);
theta=p(2);
y_trial=A*cos(x+theta);  % y with the trial parameters.
difference=y-y_trial;    % The difference between fit and experimental data
sum_square_errors = sum(difference.^2);  % The sum of the squared errors

You can put the above in a separate file called "sum_square_errors.m", or put it in your main file.

The second step is to call fmins with the appropriate arguments to minimize the sum of the squared errors. For this example I will make up some data, add noise to it and call it y. In real life, you will probably type your vectors of x and y in by hand.

x=0:.1:2*pi;
y_true = 2*cos(x+pi/4);
y_data=y_true+randn(size(x)); % Add noise to our true values
p0=[1.5 1] % Guessed parameters.
p=fmins(sum_square_errors, p0, [],[], x, y); % Use fmins to find the best fit p
plot(x,y_true,x,y_data,'r*',x,p(1)*cos(x+p(2)),'g-')
legend('True', 'Noisy', 'Best Fit')
title('Fitting Example with fmins')
xlabel('x')




function sse = expmodel(p,x,y) %EXPMODEL Fit of a 2-term exponential model to a set of data

% The parameters are in p (4-by-1 vector) and assumed to % be in the following order: a1 = p(1); b1 = p(2); a2 = p(3); b2 = p(4);

ypred = a1*exp(b1*x) - a2*exp(b2*x);  % predictions err = y - ypred;  % mismatches

sse = sum(err.^2);  % residual measure