## Weighted polynomial fitting in MATLAB (without any toolboxes)

January 7th, 2011 | Categories: matlab, programming | Tags:

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.

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.

1. I’ve been reading your blog for a while now, this fell in with something I was working on yesterday morning so I worked out the expanded math. If you are interested, it may be found here http://www.sitemodelers.com/least-squares/weighted-non-weighted-quadratic-curve-fit-to-point-set/

Cheers!

2. cool! Thanks for sharing Dave.

Mike

3. Thank you! This worked very well!

4. Thank you – this was really helpful for me. Just wanted to point out a minor typo that had me confused momentarily: where you call out your solution for the weighted quadratic under your first code box, you’re showing the unweighteed solution, instead of the myfit solution.

5. Why is this,

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