How to use Octave to Fit an Arbitrary Function with fmins

From Class Wiki
Jump to: navigation, search

Suppose we want to fit experimental data to an arbitrary function. For our purposes let that function be y = A cos(t + \theta) . We can use the fmins function from the octave package optim. We are looking for a best fit \theta and A. 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 needs to be put in a separate dot m file called model.m (the same as the function).

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

The second step is to make a script 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('model', p0, [],[], x, y_data); % 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')

Fit example.png