Say we have some simple data points and we want to fit these to
a polynomial (or a linear combination of n "regressors").
x=[1:10]'
y=1+0.5*x
plot(x,y,'.')
X=[ones(size(x)),x]
b=regress(y,X)
X*b
Of course normally there will be some error on y:
y=y+0.3*randn(size(y));
plot(x,y,'.')
b=regress(y,X)
hold on
plot(x,X*b,'r')
hold off
Let's look at a practical example:
load 050517_elnod1_edit
plot(el,v)
za=90-el;
am=1./cos(za*pi/180);
plot(am,v)
axis tight
% air mass for slab atmosphere model
z=0:80;
plot(z,1./cos(z*pi/180))
plot(am,v);
We want to write a function to regress each column versus air mass.
rgain=reg_elnod(el,v)
hist(rgain)
% normalize the gains to be equal
vn=v./(repmat(rgain,[size(v,1),1]));
plot(am,vn)
All the lines now have the same slope - equal gain.
This works fine when the model we are fitting to the data is "linear in the parameters"
but what if we want to fit a sine wave of unknown phase?
x=0:0.1:10;
p=[1,2,0.5,0.7]
y=p(1)+p(2)*sin(2*pi*(p(3)*x+p(4)));
y=y+0.1*randn(size(y));
plot(x,y,'.');
grid
Need to use iterative search algorithm.
Function lsqnonlin scans p value seeking to minimize sum(func.^2))
pf=lsqnonlin(@(v) y-sinusoid(v,x),p)
plot(x,y,'.'); hold on; plot(x,sinusoid(pf,x),'r'); hold off
Of course normally we don't know starting values and need to estimate them.
pp(1)=mean(y);
pp(2)=0.5*(max(y)-min(y));
pp(3:4)=[0.4,0];
pf=lsqnonlin(@(v) y-sinusoid(v,x),pp)
plot(x,y,'.'); hold on; plot(x,sinusoid(pf,x),'g'); hold off
This doesn't work - it wanders off to nowhere with 4 free parameters.
Can help it by fixing the parameters we know are already close to right
and doing a lower dimensional search on the others.
q=lsqnonlin(@(v) y-(pp(1)+pp(2)*sin(2*pi*(v(1)*x+v(2)))),pp(3:4));
pf=[pp(1:2),q]
plot(x,y,'.'); hold on; plot(x,sinusoid(pf,x),'g'); hold off
% refine the fit with all par free
pf=lsqnonlin(@(v) y-sinusoid(v,x),pf)
