2
$\begingroup$

I have a set of observation data {x,y} and need to make some statistical analysis. I made a density plot as follows and the blue line is the best fit of the data.

However, it seems that the blue line does not fit to the data very much. I wonder it may be affected by some noisy data (those with little counts).

Is it possible to filter the data with little counts (i.e. the data over purple area) before making the plot and best fit?

Many thank in advance!!

enter image description here

$\endgroup$
  • 4
    $\begingroup$ Don't do it! Don't get rid of "inconvenient data" nor should you smooth before fitting. A line is not necessarily the "best fit" and maybe for this data (assuming the objective is to predict the vertical axis variable from the horizontal axis variable) a more complex function is desired. Can you supply the data somewhere? $\endgroup$ – JimB Feb 10 '17 at 4:08
  • $\begingroup$ Thanks a lot. I uploaded the data to the following link. drive.google.com/file/d/0B3sHo6MIAjziaTloRnlGbm1Nam8/… $\endgroup$ – wkong Feb 10 '17 at 5:48
  • 1
    $\begingroup$ From looking at your plot and your data, I think that @AntonAntonov 's Quantile Regression package is probably your best bet. $\endgroup$ – JimB Feb 10 '17 at 6:54
  • $\begingroup$ Thank you for the suggestions. The @AntonAntonov 's Quantile Regression package looks a bit complicated to me and I need more time to digest it. Before that, I made some simple counting and filtered those with less than 20 counts on the density plot. Here is the plot. [link](drive.google.com/file/d/0B3sHo6MIAjziYVRTaER2WWVhSFE/… However, it seems that the best fit is still not the best fit. The {x,y} data is supposed to be linearly dependent in theory. I don't know what's wrong during the best fit.... $\endgroup$ – wkong Feb 10 '17 at 9:20
6
$\begingroup$

2D envelopes

One approach to get what OP wants is described in the posts

Here is code following the latter post.

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/QuantileRegression.m"]

data = Import["~/data.txt", "CSV"];

Histogram3D[data]

enter image description here

Find a threshold for what are considered outliers:

qreg = QuantileEnvelopeRegion[data, 0.98, 20];   
Show[{Graphics[{Opacity[0.02], Point[data]}, PlotRange -> All],
   BoundaryDiscretizeRegion[qreg], 
  Graphics[{Opacity[0.02], Red, Point[data]}]}, Frame -> True, 
 ImageSize -> 500]

enter image description here

Next we separate the outlier points from the rest:

rmFunc = RegionMember[qreg];

AbsoluteTiming[pred = rmFunc /@ data;]

(* {3.89161, Null} *)

Tally[pred]

(* {{True, 91123}, {False, 13330}} *)

Graphics[{Gray, Point[Pick[data, pred]], Red, 
  Point[Pick[data, Not /@ pred]]}, Frame -> True, ImageSize -> 500]

enter image description here

Note that some high count points are considered outliers. A more faithful algorithm can be devised using Quantile Regression.

Quantile regression

First we reverse a 10% sample of the data:

data2 = Map[Reverse, RandomSample[data, 10000]];

Next we select outlier thresholds and find regression quantiles:

qs = {0.1, 0.9};
AbsoluteTiming[
 qFuncs = QuantileRegression[data2, 20, qs];
 ]

(* {24.8839, Null} *)

Plot the transposed data sample and the found regression quantiles:

Show[{
  Graphics[{Opacity[0.2], Pink, Point[data2]}, PlotRange -> All],
  ListLinePlot[
   Table[{#, rq[#]} & /@ Sort[data2[[All, 1]]], {rq, qFuncs}], 
   PlotLegends -> qs]
  }, Frame -> True]

enter image description here

As in the 2D envelope procedure above, find the separation predicate:

AbsoluteTiming[
 pred2 = Map[#[[2]] >= qFuncs[[1]][#[[1]]] && #[[2]] <= 
       qFuncs[[2]][#[[1]]] &, Reverse /@ data];
]

(* {23.2092, Null} *)

Tally[pred2]

(* {{True, 83680}, {False, 20773}} *)

Plot the outlier (low counts) points and the rest:

Graphics[{Gray, Point[Pick[data, pred2]], Red, 
  Point[Pick[data, Not /@ pred2]]}, Frame -> True, ImageSize -> 500]

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ Great Thanks!!! $\endgroup$ – wkong Feb 13 '17 at 1:25
  • $\begingroup$ @wkong No problem -- good luck! $\endgroup$ – Anton Antonov Feb 13 '17 at 8:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.