Re: Optimization Problem
- From: John D'Errico <woodchips@xxxxxxxxxxxxxxxx>
- Date: Tue, 13 Sep 2005 01:24:07 GMT
In article <1126566597.791772.108540@xxxxxxxxxxxxxxxxxxxxxxxxxxxx>,
"Screwdriver223" <david.kolin@xxxxxxxxx> wrote:
> I'd like to fit data to, for lack of a better way of describing it, two
> straight lines. That is, on one interval, the data is fit by a
> straight line, and on an adjacent interval, it is fit by a second
> straight line. The "crossover" point is also a fitted parameter.
> (I've never used splines, but I think this could be described as a 1st
> order spline with 1 knot).
Almost. Implicitly you have 3 knots, even though you
may choose not to think of them that way, as you allow
the two segments to extrapolate to +/- infinity.
I prefer not to define a spline using this formulation,
but its your spline.
> I'm using lsqcurvefit with the following fitting function:
>
> function F = bilinear(x,xdata);
>
> F = zeros(size(xdata));
> F(1:round(x(5))) = xdata(1:round(x(5)))*x(1)+x(2);
> F(round(x(5))+1:end) = xdata(round(x(5))+1:end)*x(3)+x(4);
>
> x(1:4) are optimized "well", but x(5) (i.e., the crossover point)
> always returns the initial guess. I suspect this is caused by
> lsqcurvefit not taking large enough "steps" when changing x(5), but I
> don't know how to correct this. Does anyone have any suggestions?
Lets see. Two linear segments, so there are a pair of
parameters for each segment. Plus you are optimizing
the knot position. While it looks like 5 parameters,
you forget that the result will be discontinuous.
Next, by rounding the breakpoint, it will never
actually move. NEVER. You should note that it does
not do so in your code. This is why.
May I suggest a different approach? Optimize only over
the knot position, since once the knot position is
given, the spline parameters themselves are estimable
using linear least squares. Next, use a plus function
approach. (Its also NOT my preferred way to fit a
spline, but it is easy to use.)
Thus, for data in x, a breakpoint b, and (3) spline
coefficients in coef,
plusfun = @(x) x.*(x>0);
F=@(x,brk,coef) coef(1)+coef(2)*x+coef(3)*plusfun(x-brk);
First, convince yourself that F is a continuous
function across the breakpoint in brk. The function
plusfun is zero when x is negative, but returns x for
positive x. So when x is less than brk, you get
F(x) = coef(1)+coef(2)*x
i.e., a straight line. For x greater than brk, you
still get a straight line, but the line is:
F(x) = coef(1) + [coef(2)+coef(3)]*x
The only difference is the slope is different. I'll
put it all together now in Matlab:
% make some fake data.
n = 100;
x = linspace(0,1,n)';
y = sin(pi*x) + randn(size(x))/10;
% Here is the objective for lsqnonlin. I specify
% lsqnonlin on purpose.
function [res,coef] = bilinear(brk,xdata,ydata)
n = length(x);
M = [ones(n,1),xdata,(xdata-brk).*((xdata-brk)>0)];
coef = M\ydata;
res = M*coef - ydata;
% and, here is the way I'll call it from lsqnonlin:
start = mean(x);
start
xs = sort(x);
LB = xs(2);
UB = xs(end-1);
options = optimset('disp','iter','tolfun',1.e-12,'tolx',1.e-12);
brk = lsqnonlin('bilinear',start,LB,UB,options,x,y);
[junk,coef] = bilinear(brk,x,y);
% Note the trick I've used here. The linear coefficients
% were generated in bilinear. It returns two results,
% only one of which is used by lsqnonlin. The second
% one is useful in that last call to bilinear, AFTER
% the call to lsqnonlin.
Here are the spline parameters from my data:
coef =
0.0818
2.3089
-4.2668
brk =
0.4737
Do these parameters make sense? The function is
sin(pi*x) over the interval [0,1]. Its highest
curvature is at the midpoint, x = 0.5. The slopes
are quite reasonable too. Remember that to the
right of brk, the slope is the sum coef(2)+coef(3).
We can evaluate the function and plot the result
using the functions I defined above.
plusfun = @(x) x.*(x>0);
F=@(x,brk,coef) coef(1)+coef(2)*x+coef(3)*plusfun(x-brk);
plot(x,y,'b.',x,F(x,brk,coef),'r-')
HTH,
John D'Errico
--
The best material model of a cat is another, or
preferably the same, cat.
A. Rosenblueth, Philosophy of Science, 1945
.
- References:
- Optimization Problem
- From: Screwdriver223
- Optimization Problem
- Prev by Date: Re: Visibility of ActiveX Control in GUI
- Next by Date: Re: Install problem: error while loading shared li
- Previous by thread: Optimization Problem
- Next by thread: Re: Optimization Problem
- Index(es):
Relevant Pages
|
Loading
