# Finding local extrema in noisy data using Quantile Regression

## Introduction

This blog post (article) describes an algorithm for finding local extrema in noisy data using Quantile Regression. The problem formulation and a solution for it using polynomial model fitting (through LinearModelFit) were taken from Mathematica StackExchange — see “Finding Local Minima / Maxima in Noisy Data”, [1].

The proposed Quantile Regression algorithm is a version of the polynomial fitting solution (proposed by Leonid Shifrin in [1]) and has the following advantages: (i) requires less parameter tweaking, and (ii) more importantly it is more robust with very noisy and oscillating data. That robustness is achieved by using two regression fitted curves: one close to the local minima and another close to the local maxima, computed for low and high quantiles respectively. (Quantile Regression is uniquely able to do that.)

The code for this blog post is available at [6].

A complete version of this blog post is available as a PDF document here: Finding local extrema in noisy data using Quantile Regression.pdf.

## The problem

Here is the problem formulation from [1].

Problem 1: We have a list of pairs of numbers representing measurements of a scalar variable on a regular time grid. The measurements have noise in them.

As a data example for this problem the author of the question in [1] has provided the following code:

```temptimelist = Range[200]/10;
tempvaluelist = Sinc[#] &@temptimelist + RandomReal[{-1, 1}, 200]*0.02;
data1 = Transpose[{temptimelist, tempvaluelist}];
ListPlot[data1, PlotRange -> All, Frame -> True, ImageSize -> 500]```

In this article we are going to consider a more general problem hinted in the discussion of the solutions in [1].

Problem 2: Assume that the data in Problem 1 is collected several times and the noise is present in both the measured values and the time of measurement. Also the data can be highly oscillatory in nature.

Consider this data generation as an example for Problem 2:

```n = 1000;
xs = N@Rescale[Range[n], {1, n}, {0, 60}];
data2 = Flatten[
Table[Transpose[{xs + 0.1 RandomReal[{-1, 1}, Length[xs]],
Map[5 Sinc[#] + Sin[#] + 4 Sin[1/4 #] &, xs] +
1.2 RandomVariate[SkewNormalDistribution[0, 1, 0.9],
Length[xs]]}], {10}], 1];
ListPlot[data2, PlotRange -> All, Frame -> True, ImageSize -> 500]```

Note that this data has 10000 points and it is much larger than the data for Problem 1.

## Extrema location approximation by model fitting followed by nearest neighbors search

Several solutions are given in [1]. A couple are using wavelets or Gaussian filtering for de-noising. The one we are going to focus on in this article is using polynomial fitting for extrema location approximation and then finding the actual data extrema by Nearest Neighbors (NN’s) search. It is provided by Leonid Shifrin.

Let us list the steps of that algorithm:

1. Fit a polynomial through the data (using LinearModelFit).
2. Find the local extrema of the fitted polynomial. (We will call them fit estimated extrema.)
3. Around each of the fit estimated extrema find the most extreme point in
the data by nearest neighbors search (by using Nearest).

We are going to refer to this algorithm as LMFFindExtrema. Its implementation is available here: QuantileRegressionForLocalExtrema.m, [6].

Here is the application of the algorithm to the data example of Problem 1:

(The continuous line shows the fitted curve.)

Here is the application of the algorithm to the data example of Problem 2:

It does not help to just increase the number of basis polynomials and the number of NN’s examined points:

We can also see that increasing the number of examined NN’s for each of the fit estimated extrema would make some the final result points to be “borrowed” from neighboring peaks in the data.

Experiments with fitting trigonometric basis functions showed very good fit to the data but the calculations of the fit estimated extrema were very slow.

## Using Quantile Regression

It is trivial to rewrite the algorithm LFMFindExtrema to use Quantile Regression instead of linear model fitting of polynomials. We are going to use the Mathematica package [2] implementation described in [3,4]. The function QuantileRegression provided by [2] uses a B-spline basis [5] to find the quantile regression curves (also known as regression quantiles).

We are going to call this algorithm QRFindExtrema. Again its implementation is available here: QuantileRegressionForLocalExtrema.m, [6].

The algorithm QRFindExtrema has the following parameters: the data, number of B-spline knots, interpolation order, quantiles (corresponding to the curves to be fitted).

QRFindExtrema returns a list of regression quantile functions and a list of lists with extrema estimates.

More importantly though, QRFindExtrema can use two curves for finding the local extrema: one for local minima, and one for local maxima. (This feature is justified below.)

Here is the application of QRFindExtrema to the example data of Problem 1:

Here is application of QRFindExtrema to the data of Problem 2:

We can see that the regression quantile for 0.5 is too flat to get good judgment of the local extrema. We can get better results if we increase the number of knots for the B-spline basis built by QuantileRegression.

We can see that results look “almost right”, the horizontal locations of the peaks are apparently more-or-less correctly identified, but the result extrema are too close to the fitted curve. Just increasing the number of examined NN’s for the fit estimated extrema does not produce good results because points from neighboring peaks are being chosen as final extrema estimates.

In order to solve this problem we use two regression quantiles. For local minima we use a regression quantile for a low quantile number, say, 0.02; for the local maxima we use a regression quantile for a large quantile number, say, 0.98 .

The use of two (or more) curves to be fitted is an unique capability of Quantile Regression. With this algorithm feature by construction the lower regression quantile is close to the local minima and the higher regression quantile is close to the local maxima.

Also, since we find two regression quantile curves we can use two nearest neighbors finding functions: one with the points below the low regression quantile, and one with the points above the high regression quantile. The implementation in [6] takes an option specification for should the nearest neighbor functions for finding the extrema be constructed using all data points or just the outliers (the points outside of the found regression quantiles).

## More examples

More timid noisy and oscillating data with around 10,000 points.

## References

[1] Mathematica StackExchange discussion. “Finding Local Minima / Maxima in Noisy Data”,
URL: http://mathematica.stackexchange.com/questions/23828/finding-local-minima-maxima-in-noisy-data/ .

[2] Anton Antonov, Quantile regression Mathematica package, source code at GitHub, https://github.com/antononcube/MathematicaForPrediction, package QuantileRegression.m, (2013).

[3] Anton Antonov, Quantile regression through linear programming, usage guide at GitHub, https://github.com/antononcube/MathematicaForPrediction, in the directory “Documentation”, (2013).

[4] Anton Antonov, Quantile regression through linear programming, “Mathematica for prediction algorithms” blog at WordPress.com, 12/16/2013.

[5] Anton Antonov, Quantile regression with B-splines, “Mathematica for prediction algorithms” blog at WordPress.com, 1/1/2014.

[6] Anton Antonov, QuantileRegressionForLocalExtrema Mathematica package, source code at GitHub, https://github.com/antononcube/MathematicaForPrediction, package QuantileRegressionForLocalExtrema.m, (2015).
The package is in the directory “Applications”.

# Directional quantile envelopes in 3D

## Introduction

This blog post was mostly made as a response to a comment of my previous blog post “Directional quantile envelopes”, [1]:

This looks extremely useful and elegant – any plans on generalizing to 3 (or more) dimensions?
–Jan Liphardt

Since I did say in the previous blog post that the algorithm can be easily extended to 3D and higher dimensions, I decided to illustrate that with an example. The 3D implementation is really easy using Mathematica’s geometric computation functions (i) for derivation of rotation matrices and (ii) for regions manipulation. (The latter are new in version 10.)

The algorithm shown below is a slightly modified version of the algorithm in “Directional quantile envelopes”. Instead of the step that does the intersection of directional quantile lines with Reduce, we can use ImplicitRegion to specify the region enclosed by the directional quantile planes. We can use other built-in functions for region predicates and manipulation to derive, visualize, or utilize the directional quantile regions found by the algorithm. (All these are exemplified below.)

## Data generation

Let us first generate some data.

```npoints = 10000; data = {RandomReal[{0, 2 Pi}, npoints], RandomReal[{-Pi, Pi}, npoints], RandomVariate[PoissonDistribution[4.8], npoints]}; data = MapThread[#3*{Cos[#1] Cos[#2], Sin[#1] Cos[#2], Sin[#2]} &, data];```

Let us plot the data

```Block[{qs = 12}, qs = Map[Quantile[#, Range[0, 1, 1/(qs - 1)]] &, Transpose[data]]; ListPointPlot3D[data, PlotRange -> All, PlotTheme -> "Detailed", FaceGrids -> {{{0, 0, -1}, Most[qs]}, {{0, 1, 0}, qs[[{1, 3}]]}, {{-1, 0, 0}, Rest[qs]}}, ImageSize -> 700] ]```

On the plot the grid lines are derived from the quantiles of x, y and z coordinates of the data set.

## Algorithm application step-by-step

1. Standardize the data
1.1. This step is generally not needed and in this implementation would complicate things.

2. Create a set of uniform angles. (Here we generate vectors that define 3D directions.)

```nqs = 20; dirs = N@Flatten[ Table[{Cos[\[Theta]] Cos[\[Phi]], Sin[\[Theta]] Cos[\[Phi]], Sin[\[Phi]]}, {\[Theta], 2 \[Pi]/(10 nqs), 2 \[Pi], 2 \[Pi]/nqs}, {\[Phi], -\[Pi], \[Pi], 2 \[Pi]/nqs}], 1];```

`Graphics3D[Map[Arrow[{{0, 0, 0}, #}] &, dirs]]`

3. Rotate the data for each direction vector, a ∈ A. (Here we just generate the rotation matrices.)

```In[672]:= rmats = RotationMatrix[{{1, 0, 0}, #}] & /@ dirs; Length[rmats]```

``` ```

`Out[673]= 420`

4. Find the q-quantile of the z-coordinates of the data rotated at direction a. Denote that point with Z[a,q].

```In[674]:= qs = {0.95}; qDirPoints = Flatten[Map[Function[{m}, Quantile[(m.Transpose[data])[[3]], qs]], rmats]];```

5. Find the points z[a,q] corresponding to Z[a,q] in the original coordinate system.
5.1. In this implementation this step is redundant, see the next step.

6. Find the quantile planes intersections
6.1. Instead of doing this step we simply use ImplicitRegion.

```In[676]:= qRegion = ImplicitRegion[ MapThread[(#1.{x, y, z})[[3]] <= #2 &, {rmats, qDirPoints}], {x, y, z}];```

`In[677]:= Shallow[qRegion]`

``` ```

```Out[677]//Shallow= ImplicitRegion[ LessEqual[<>] && LessEqual[<>] && LessEqual[<>] && LessEqual[<>] && LessEqual[<>] && LessEqual[<>] && LessEqual[<>] && LessEqual[<>] && LessEqual[<>] && LessEqual[<>] && <>, {x, y, z}]```

## Wrapping it in a function

Here is a definition of a function that wraps all of the algorithms steps in the previous section.

``` Clear[QuantileEnvelopeRegion] QuantileEnvelopeRegion[points_?MatrixQ, quantile_?NumberQ, numberOfDirections_Integer] := Block[{nd = numberOfDirections, dirs, rmats, qDirPoints, qRegion}, dirs = N@Flatten[ Table[{Cos[\[Theta]] Cos[\[Phi]], Sin[\[Theta]] Cos[\[Phi]], Sin[\[Phi]]}, {\[Theta], 2 \[Pi]/(10 nd), 2 \[Pi], 2 \[Pi]/nd}, {\[Phi], -\[Pi], \[Pi], 2 \[Pi]/nd}], 1]; rmats = RotationMatrix[{{1, 0, 0}, #}] & /@ dirs; qDirPoints = Flatten[Map[ Function[{m}, Quantile[(m.Transpose[points])[[3]], quantile]], rmats]]; qRegion = ImplicitRegion[ MapThread[(#1.{x, y, z})[[3]] <= #2 &, {rmats, qDirPoints}], {x, y, z}]; qRegion ] /; Dimensions[points][[2]] == 3 && 0 < quantile <= 1; ```

## Visualizing

From the implicit region specification we can generate an approximation of the boundary surface using the built-in function BoundaryDiscretizeRegion.

```qRegion = QuantileEnvelopeRegion[data, 0.95, 20]; qRegion2 = QuantileEnvelopeRegion[data, 0.8, 20];```

```qSurface = BoundaryDiscretizeRegion[qRegion]; qSurface2 = BoundaryDiscretizeRegion[qRegion2]; ```

`Grid[{{qSurface, qSurface2}}]`

Now we can visualize the quantile surface together with the data points:

Block[{c = 3, opts, pntsgr},
``` opts = {PlotRange -> {{-c, c}, {0, c}, {-c, c}}, Boxed -> True, Axes -> True, ImageSize -> {Automatic, 500}}; pntsgr = Graphics3D[Point[data]]; Grid[{{Show[{qSurface, pntsgr}, opts], Show[{qSurface2, pntsgr}, opts]}}] ]```

## Check

Now we can answer the question how many data points are inside the quantile directions regions. Again, this is easy with the built-in region functionalities of Mathematica version 10.

```In[656]:= inPoints = Select[data, # \[Element] qRegion &]; Length[inPoints]/Length[data] // N```

``` ```

`Out[657]= 0.5974`

```In[700]:= inPoints = Select[data, # \[Element] qRegion2 &]; Length[inPoints]/Length[data] // N```

``` ```

`Out[701]= 0.1705`

(In these code lines I kept “[\Element]” instead of replacing it with “∈” in order the lines to be copy and pasted. )
As we can see the algorithm would produce a region Subscript[R, q] which contains a much smaller number of points than the requested quantile, q. This property of the algorithm is discussed in [2]. (In the two dimensional space this property is less pronounced.)

Obviously, though, we would expect the bias of the algorithm to be monotonic with respect to the quantiles requested. Given a set of points P and quantiles u and v, 0<u<v<=1, for the respectively produced by the algorithm regions Ru and Rv we have |P ∩ Ru| < |P ∩ Rv| .

## References

[1] Anton Antonov, Directional quantile envelopes, “Mathematica for prediction algorithms” blog at WordPress.com, 11/3/2014 .

[2] Linglong Kong, Ivan Mizera, “Quantile tomography: using quantiles with multivariate data”, 2013, arXiv:0805.0056v2 [stat.ME] URL: http://arxiv.org/abs/0805.0056 .

# Quantile regression with B-splines

Last weekend I made two implementations of Quantile Regression (QR) calculation with B-spline bases. The first implementation is based on the Linear Programming (LP) formulation of the quantile minimization problem. The second implementation is a direct translation of the non-LP minimization formulation. (Mathematica’s functions LinearProgramming and Minimize are used respectively.)

The implementations are in the package QuantileRegression.m provided by the project MathematicaForPrediction at GitHub.

Because of the signature required by the QR implementation with B-splines and the returned results, I had to rename the previous QR implementation to QuantileRegressionFit and I named the new one (with the B-splines) QuantileRegression. I think the new names reflect better what the functions are doing.

Right now, to me the most interesting application of QR is detection of outliers in time series and reconstruction of conditional density functions. Below is a demonstration of using QR for outlier detection.

Consider this data:

which was generated with the code:
``` tempData = WeatherData["Atlanta", "Temperature", {{2006, 1, 1}, {2013, 12, 31}, "Day"}]; tempPData = tempData; tempPData[[All, 1]] = DateDifference[tempData[[1, 1]], #, "Day"] & /@ tempData[[All, 1]]; tempPData[[All, 1]] = tempPData[[All, 1, 1]]; tempPData[[All, 2]] = Sqrt[tempPData[[All, 1]]] + tempPData[[All, 2]]; ```
(I took the temperature in Atlanta, GA, USA for the last 7 years and modified the data with the formula Yi = sqrt(Xi) + Yi .)

Here is a plot with five regression quantiles calculated using B-splines for the quantiles {0.02, 0.2, 0.4, 0.6, 0.8, 0.98}:

Here is another plot with the 0.02 and 0.98 regression quantiles outlining the data:

The regression quantiles were computed with the command:
``` qFuncs = QuantileRegression[tempPData, 44, {0.02, 0.2, 0.4, 0.6, 0.8, 0.98}] ```

Note that in order to get satisfactory results, the only tuning value I had to choose is for the second parameter, which in this case specifies the number of B-spline basis knots. The second parameter can also be a list of knots. The option InterpolationOrder specifies the splines order. (By default InterpolationOrder->3.)

### Update, 2015.07.04

This blog post was written with an older version of Mathematica, version 9. (The current is 10.1.)  Mathematica 10 uses different format for the `WeatherData` results — `TimeSeries` object are returned, which use `QuantityArray` and/or `Quantity` around the numerical values like temperature, wind speed, etc.

Below is given code that runs with Mathematica 10.1. (The only change is the second line of the data loading part.)

``` tempData = WeatherData["Atlanta", "Temperature", {{2006, 1, 1}, {2013, 12, 31}, "Day"}]; tempData = Transpose[{Normal[tempData["Dates"]], Normal[tempData["Values"]] /. Quantity[q_, _] :> q}]; tempPData = tempData; tempPData[[All, 1]] = DateDifference[tempData[[1, 1]], #, "Day"] & /@ tempData[[All, 1]]; tempPData[[All, 1]] = tempPData[[All, 1, 1]]; tempPData[[All, 2]] = Sqrt[tempPData[[All, 1]]] + tempPData[[All, 2]]; ```

2. Regression quantiles computation:
``` qFuncs = QuantileRegression[tempPData, 44, {0.02, 0.2, 0.4, 0.6, 0.8, 0.98}]; ```

3. Plots:
``` gr1 = ListPlot[tempPData]; gr2 = Plot[ Evaluate[Through[qFuncs[x]]], {x, Min[tempPData[[All, 1]]], Max[tempPData[[All, 1]]]}, PlotPoints -> 1200, PerformanceGoal -> "Speed"]; jan1Pos = Position[DateList /@ tempData[[All, 1]], {_, 1, 1, ___}][[All, 1]]; Show[{gr1, gr2}, PlotRange -> All, Axes -> False, Frame -> True , GridLines -> {jan1Pos, FindDivisions[{##}, 10] &}, GridLinesStyle -> Directive[GrayLevel[0.8], Dashed], AspectRatio -> 1/3] ```

# Quantile regression through linear programming

We can say that least squares linear regression corresponds to finding the mean of a single distribution. Similarly, quantile regression corresponds to finding quantiles of a single distribution. With quantile regression we obtain curves — “regression quantiles” — that together with the least squares regression curve would give a more complete picture of the distributions (the Y’s) corresponding to a set of X’s.

The Quantile Regression Problem (QRP) is formulated as a minimization of the sum of absolute differences. Further, QRP is re-formulated as a linear programming problem, which allows for efficient computation. To compute quantiles other than the median the so called “tilted” function is used. (See this Wikipedia entry.)

For a complete, interesting, and colorful introduction and justification to quantile regression I refer to this article:

Roger Koenker, Gilbert Bassett Jr., “Regression Quantiles”, Econometrica, Vol. 46, No. 1. (Jan., 1978), pp. 33-50. JSTOR URL: http://links.jstor.org/sici?sici=0012-9682%28197801%2946%3A1%3C33%3ARQ%3E2.0.CO%3B2-J .

I recently implemented and uploaded a Mathematica package for computing regression quantiles using Mathematica’s function LinearProgramming — see the package QuantileRegression.m provided by the MathematicaForPrediction project at GitHub. Also see this guide for using the function QuantileRegressionFit provided by the package. The guide has some theoretical explanations and shows how the quantile regression problem can be formulated as a linear programming problem.

Below are presented several examples of regression quantiles computed over different data sets with the function QuantileRegressionFit provided by the package. The signature of QuantileRegressionFit is very similar to that of Mathematica‘s function Fit. Here is the definition statement:

QuantileRegressionFit::usage = “QuantileRegressionFit[data,funs,var,qs] finds the regression quantiles corresponding to the quantiles qs for a list of data as linear combinations of the functions funs of the variable var.”

QuantileRegressionFit has a Method option and can use both Minimize and LinearProgramming, but the computations with Minimize are quite slow for larger data sets. (The implementation with Minimize is included for didactic purposes.) All QuantileRegressionFit results presented below are done with the linear programming implementation.

First let us consider data generated by adding skewed noise over a logarithmic curve:

Pretending that we don’t know how the data is generated, just by looking at the plot, we might conjecture that the model for the data is

Y = b0 + b1 X + b2 X^(1/2) + b3 log(X) .

Using this model we find regression quantiles corresponding to 0.05, 0.25, 0.5, 0.75, 0.95. For comparison we also find the least squares fit of the model. Here is a plot that shows the regression quantiles for that data together with the least squares fit of the model functions:

We can check do the regression quantiles partition the data in the expected way. Here is a table that shows the fraction of the data points with Y’s greater or equal than the corresponding regression quantile:

Here is a plot of the with the regression quantiles of a more interesting set of data:

Here is the table for the data partition tests:

I further made some performance profiling tests. First we need to choose a family or several families of test data. Also, since Mathematica’s function LinearProgramming has several methods it is a good idea to test with all of them. Here I am going to show results only with one family of data and two LinearProgramming methods. The data family is the skewed noise over a logarithmic curve used as an example above. The first LinearProgramming method is Mathematica‘s  (default) “InteriorPoint”, the second method is “CLP” that uses the built-in COIN-OR CLP optimizer. I run the profiling tests using one quantile {0.5} and five quantiles {0.05, 0.25, 0.5, 0.75, 0.95}, which are shown in blue and red respectively. I also run tests with different number of model functions {1, x, x^(1/2), log(x)} and {1, x, log(x)} but there was no significant difference in the timings (less than 2%).

It is interesting to note that the average ratio of the timings with 1 vs. 5 quantiles is 0.38 for “InteriorPoint” and 0.5 for “CLP”.

In another post I am going to show examples of the robustness of regression quantiles when the data has outliers.

Update: Because I implemented calculation of regression quantiles with B-splines I had to rename the function QuantileRegression into QuantileRegressionFit. I made the corresponding changes in the blog post. (Except the profile graphics captions.)