Basically I have a two dimensional list plot with a bunch of points between the square of (0,0) and (1,1). The y-values of these points are drawn from a uniform probability distribution between 0 and 1, and the x-values are drawn from a semi-uniform distribution, except it jumps up at x=.5. So in the plot the density of points jumps up at x=0.5. I need a statistical technique that will determine where this jump occurs and will draw a line there, even if we didn't know the jump was at 0.5 or if the initial conditions were changed and it wasn't at 0.5 or even a straight line.

List Plot

  • 1
    $\begingroup$ Off the top of my head, subdivide the horizontal axis in small steps (bins), count the number of points that fall into each bin (let's call that 'density') and find where the value of density has a jump higher than your selected threshold. I believe there's a function named Binning or Bincount to do most of the work. $\endgroup$ – Peltio Nov 20 '13 at 3:50
  • $\begingroup$ What code have you already tried? $\endgroup$ – Hector Nov 20 '13 at 4:25
  • $\begingroup$ You can make a SmoothKernelDistribution, plot its PDF, visually find the approximate density (PDF value) at the jump, then use FindRoot to find the location of the jump. $\endgroup$ – Szabolcs Nov 20 '13 at 4:56

I think there are standard statistical methods -- there are statisticians that visit this site who might help. Otherwise, you might ask on the stats.SE site for an algorithm.

Here's a naive way: Find the jump in the slope of the CDF.

data = RandomReal[1, {1000, 2}] ~Join~
   Transpose[{0.5, 0} + {0.5, 1.} RandomReal[1, {2, 1000}]];
distx = EmpiricalDistribution[data[[All, 1]]];

Plot[CDF[distx, x], {x, 0, 1 - 0.02}, PlotRange -> All]

Mathematica graphics

The second differences show a peak at the jump. You may have to tune the dx parameter by hand, but one might be able to estimate it from statistics of the data set. It seems to produce reasonable estimates of where the increase in the density occurs for a fairly wide range of values of dx. The estimate varies by about 5%+ (of 0.5).

With[{dx = 0.03},
 Plot[Differences[CDF[distx, x + {-dx, 0, dx}], 2], {x, 0, 1 - dx}, PlotRange -> All]

Mathematica graphics

Find the maximum:

Block[{f, dx = 0.03},
 f[x_?NumericQ] := 
  First@Differences[CDF[distx, x + {-dx, 0, dx}], 2];
 NMaximize[{f[x], 0 <= x <= 1}, x]

(* {0.0355, {x -> 0.502374}} *)


ListPlot[data, GridLines -> 
   {{x} /. Last[Block[{f, dx = 0.03}, 
       f[(x_)?NumericQ] := First[Differences[
           CDF[distx, x + {-dx, 0, dx}], 2]]; 
        NMaximize[{f[x], 0 <= x <= 1}, x]]], None}, 
  GridLinesStyle -> Red]

Mathematica graphics

What would be nice is a method for computing the confidence interval (either of this method or another). Unfortunately that is beyond my understanding of statistics.

| improve this answer | |
  • $\begingroup$ Thank you!! I think that's exactly what I needed. I might be back to ask how to do it for a line that isn't necessarily vertical. $\endgroup$ – Ryan Summers Nov 21 '13 at 4:11

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.