MATLAB Central Newsreader - An Interesting Real-World Problem...

An Interesting Real-World Problem...

Geoffrey Akien — Thu, 02 Jul 2009 18:34:01 -0400

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

Re: An Interesting Real-World Problem...

Adam Chapman — Thu, 02 Jul 2009 18:54:16 -0400

On Jul 2, 7:34 pm, "Geoffrey Akien" <geoff.ak...@gmail.com> wrote:
> 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

is the sharp drop always at the same gradient? Im thinking a step
function might be adequate

Re: An Interesting Real-World Problem...

Adam Chapman — Thu, 02 Jul 2009 18:55:40 -0400

On Jul 2, 7:34 pm, "Geoffrey Akien" <geoff.ak...@gmail.com> wrote:
> 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

Is the gradient of the sharp drop always the same? Im thinking a
macauley step function might be useful

Re: An Interesting Real-World Problem...

Stan Bischof — Thu, 02 Jul 2009 19:00:04 -0400

Geoffrey Akien <geoff.akien@gmail.com> wrote:
> 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)
>

If you are looking for the "edge" and your data is always about the same shape
as the picture that you have then I would suggest that your best answer
would be the same as you would get from pencil-and-paper:

I would draw three lines to model this data on paper, and the intersections
of those three lines would bound the regions. The three lines would be the
best fits in the three regions so in this case a upward sloped line,
then a downward sloped line and then a flat line.

Do exactly the same on your computer by modeling it with a
piece-wise linear curve fit. You could chunk out the least-squares
analytic solution but there's likely a piecewise linear
function in Matlab curve-fitting toolbox.

A linear spline fit would be about the same thing and then your control
points would be the corners.

Stan

Re: An Interesting Real-World Problem...

Rune Allnor — Thu, 02 Jul 2009 19:16:59 -0400

On 2 Jul, 20:34, "Geoffrey Akien" <geoff.ak...@gmail.com> wrote:
> 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?

Don't know about 'good', but at least a suggestion:

1) Use a 5-7 pt median filter to get rid of noise
2) Differentiate result
3) Integrate twice
4) Determine the threshold as where the cumulative
sum reaches, say, 99% of max.

Or something like that.

Rune

Re: An Interesting Real-World Problem...

TideMan — Thu, 02 Jul 2009 20:11:49 -0400

On Jul 3, 6:34 am, "Geoffrey Akien" <geoff.ak...@gmail.com> wrote:
> 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 orthogonal wavelet decomposition. Empirical mode
decomposition is another alternative.
Here's a web page where I show how you can apply wavelet decomposition
to a signal similar to yours to extract information about the signal
at different timescales:
http://www.tideman.co.nz/Salalah/OrthWaveDecomp.html

Re: An Interesting Real-World Problem...

John D'Errico — Thu, 02 Jul 2009 21:09:02 -0400

"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

Re: An Interesting Real-World Problem...

Geoffrey Akien — Fri, 03 Jul 2009 17:25:03 -0400

Alright! SLMtools worked a treat with the piece-wise linear fit - I got exactly what I wanted out of it. Thanks!

Re: An Interesting Real-World Problem...

John D'Errico — Fri, 03 Jul 2009 19:17:01 -0400

"Geoffrey Akien" <geoff.akien@gmail.com> wrote in message <h2letf$nbb$1@fred.mathworks.com>...
> Alright! SLMtools worked a treat with the piece-wise linear fit - I got exactly what I wanted out of it. Thanks!

Good. If you find you need to provide starting values
for the knots to get convergence, just pass them in as
a vector for the knots parameter.

John