Path: news.mathworks.com!not-for-mail
From: "John D'Errico" <woodchips@rochester.rr.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: An Interesting Real-World Problem...
Date: Thu, 2 Jul 2009 21:09:02 +0000 (UTC)
Organization: John D'Errico (1-3LEW5R)
Lines: 81
Message-ID: <h2j7le$jmi$1@fred.mathworks.com>
References: <h2iuip$nm2$1@fred.mathworks.com>
Reply-To: "John D'Errico" <woodchips@rochester.rr.com>
NNTP-Posting-Host: webapp-03-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1246568942 20178 172.30.248.38 (2 Jul 2009 21:09:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Thu, 2 Jul 2009 21:09:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 869215
Xref: news.mathworks.com comp.soft-sys.matlab:552554


"Geoffrey Akien" <geoff.akien@gmail.com> wrote in message <h2iuip$nm2$1@fred.mathworks.com>...
> I've been using MATLAB for around 6 months now, and only now have I needed to do something that is actually mathematically-related!  I'm a chemist, so high-level maths is not my strong point.
> 
> I have some real-world data that I'd like to pick out the dramatic change in gradient of:
> 
> http://www.flickr.com/photos/40058045@N07/3681613485/
> (wish I could include this in the post but never mind)
> 
> Because the data is quite noisy (and not random) this is quite awkward.  In that picture there are about 5000 points.  You cannot smooth too much (in practice >300) without losing the information about the point at which the gradient changes (the smoothed data starts to cut the corner).  Even when the data is smoothed, the presence of non-random noise means that the "flat" part of the curve is not actually flat, and there are several stationary points (calculated thanks to STATIONARY from FEX) along that flat part.
> 
> I've also tried fitting a least-squares linear regression using POLYFIT for each of the data points with ~500 points either side of it.  This was with the aim of picking out the part of the graph with the lowest R2 value - this I assumed would represent the fast gradient change.  However, this still does not pick out the right point, selecting a point about 0.1 hours before it should be, and I can't think why this is the case.
> 
> Any other good suggestions?
> 
> Thanks
> 
> Geoff

I'd use my slmtools. It is a least squares spline
tool that can be constrained to have various
properties.

You can force it to be a constant function above
a certain point, and monotone decreasing above
another point. Throw enough knots at the curve
to be able to resolve the break, and then constrain
it to behave as you know it should.

Alternatively, use the ability of this tool to fit a
piecewise linear least squares function, where it
will then choose the location of the break points
to best fit the data. It can use all of the abilities
to constrain the curve shape as the cubic version.

The call in the latter case might look like 

model = slmengine(x,y,'knots',4,'rightslope',0, ...
      'degree',1,'interiorknots','free','plot','on');

1. This fits the curve as a piecewise linear spline,
with exactly 3 segments (therefore 4 knots).

2. It forces the piece on the right end to have a
zero slope, i.e. it will be a constant segment.

3. The placement of the middle two knots will be
chosen by the tool itself. (It uses fmincon, so you
will need the optimization toolbox.)

4. It plots the data and the curve overlaid upon it.

To fit the curve using a cubic spline, I might have
tried something vaguely like this instead,

model = slmengine(x,y,'knots',20,'rightslope',0, ...
      'degree',3,'plot','on','decreasing',[1.3,max(x)], ...
      'concavedown',[min(x),1.4]);

1. This fits a curve that has 20 knots, equally spaced.
Just use a lot of knots here, then as I said, constrain
the curve to do as you expect.

2. Build it as a C2 cubic spline.

3. Plot the results.

4. Make the curve a monotone decreasing function
above x = 1.3.

5. Force it to have a non-positive second derivative
below x = 1.4.

As I said, it will require the optimization toolbox.
Find SLM for download from the file exchange. There
are an extensive set of examples to show how it can
be used there.

http://www.mathworks.com/matlabcentral/fileexchange/24443

HTH,
John