Walking Randomly

Someone recently contacted me with a problem – she wanted to use MATLAB to perform a weighted quadratic curve fit to a set of data.  Now, if she had the curve fitting toolbox this would be nice and easy:

x=[ 1 2 3 4 5];
y=[1.4 2.2 3.5 4.9 2.3];
w=[1 1 1 1 0.1];

f=fittype('poly2');
options=fitoptions('poly2');
options.Weights=[1 1 1 1 0.1];
fun=fit(x',y',f,options);
p1=fun.p1;
p2=fun.p2;
p3=fun.p3;

myfit = [p1 p2 p3]

myfit =
   -0.1599    1.8554   -0.5014

So, our weighted quadratic curve fit is y = -0.4786*x^2 + 3.3214*x – 1.84

So far so good but she didn’t have access to the curve fitting toolbox so what to do? One function that almost meets her needs is the standard MATLAB function polyfit which can do everything apart from the weighted part.

polyfit(x,y,2)
ans =
   -0.4786    3.3214   -1.8400

which would agree with the curve fitting toolbox if we set the weights to all ones. Sadly, however, we cannot supply the weights to the polyfit function as it currently stands (as of 2010b). My solution was to create a new function, wpolyfit, that does accept a vector of weights:

wpolyfit(x,y,2,w)
ans =
   -0.1599    1.8554   -0.5014

So where is this new function I hear you ask? Normally this is where I would provide you with a download link but unfortunately I created wpolyfit by making a very small modification to the original built-in polyfit function and so I might be on dicey ground by distributing it. So, instead I will give you instructions on how to make it yourself.  I did this using MATLAB 2010b but it should work with other versions assuming that the polyfit function hasn’t changed much.

First, open up the polyfit function in the MATLAB editor

edit polyfit

change the first line so that it reads

function [p,S,mu] = wpolyfit(x,y,n,w)

Now, just before the line that reads

% Solve least squares problem.

add the following

w=sqrt(w(:));
y=y.*w;
for j=1:n+1
    V(:,j) = w.*V(:,j);
end

Save the file as wpolyfit.m and you are done. I won’t go over the theory of how this works because it is covered very nicely here.  Comments welcomed.