How can one efficiently find the (x,y) values of the rectangles that enclose a chosen % of nonrandom data points in mathematica?

Question

A half hour of searching for answers did not give useful results. This is typical data:

I need to find the (x,y) values of at least 6 rectangles with one corner at (0,0) and the opposite corner at (x,y) which enclose a defined percentage of all 10^6 data points.

4/20/18 Update: I have added Carl Woll's code to the program. It gives a fast and accurate answer for the (x,y) values of the rectangles that enclose a chosen % of the data points. Here is the full code:

(*----------------------------------------------------------------*)
(*Remove all symbols*)
Remove["Global`*"];
(*--------------------------------*)
(*Side force data in kgF*)
forceside = 
  Round[{1.46372184`, 1.5154776`, 1.62452304`, 1.64008152`, 
    1.64888136`, 1.64924424`, 1.7485826400000002`, 
    1.9107900000000002`, 1.92067848`, 1.95651288`, 
    1.9618200000000001`, 2.02047048`, 2.0909145600000003`, 
    2.13232824`, 2.1611772`, 2.22912648`, 2.36516112`, 2.43143208`, 
    2.5262344800000003`, 2.57454288`, 2.61110304`, 2.66857416`, 
    2.7966708`, 2.93379408`, 3.24033696`}, 0.0001];
(*Assume rear force data is identical to side force data*)
forcerear = forceside;
(*Bottom force data in kgF*)
forcebottom = 
  Round[{1.67655096`, 1.73234376`, 1.8965016000000001`, 2.06723664`, 
    2.19692088`, 2.20885056`, 2.21420304`, 2.2199184`, 2.25389304`, 
    2.2695422400000003`, 2.42848368`, 2.62902024`, 2.65727952`, 
    2.76596208`, 2.95570296`, 3.15574056`, 3.25290168`, 
    3.5498735999999997`, 3.8103760799999997`, 3.88707984`, 
    4.00737456`, 4.1241312`, 4.34081592`, 4.56303456`, 4.57950024`}, 
   0.0001];
(*--------------------------------*)
(*set the Monte Carlo sample size *)
mccount = 10^6;
(*--------------------------------*)
(*Sample the pull locations*)
pulllocations = RandomChoice[{0.70, 0.16, 0.14} -> {0, 0, 1}, mccount];
(*Remove all the pull locations where the forces = 0*) 
nonzeropulllocations = DeleteCases[pulllocations, 0];
(*Sample the HDMI port orientations*) 
portorientations = 
  RandomChoice[{0.41206, 0.46231, 0.12563} -> IdentityMatrix[3], 
   Length[nonzeropulllocations]];
(*Separate out three lists for the three port orientations with \
nonzero pull forces and then remove the zeroes corresponding to pulls \
on other ports*)
portside = DeleteCases[portorientations[[All, 1]], 0];
portrear = DeleteCases[portorientations[[All, 2]], 0];
portbottom = DeleteCases[portorientations[[All, 3]], 0];
(*Sample the total number of cycles over 5 years for each customer \
with nonzero pulls*)
(*The five year bins are:*) 
fiveyearbins = {{0, 0}, {1, 109}, {110, 209}, {210, 409}, {410, 
    609}, {610, 809}, {810, 1009}, {1010, 1400}};
(*Assign the relative probabilities for a customer to be in each bin from the survey data*)
userprobabilities = {43, 36, 3, 1, 1, 1, 2, 1};
(*set up eight bins (1,2,3,4,5,6,7,8} for one user to randomly fall \
into*)
oneuser = IdentityMatrix[8];
(*Sample the number of users falling into each five year bin*)
cyclenumbers =
  ParallelTable[
   (*Sample the numbers of users in each bin, for 88 users:*) 
   usersperbin = RandomChoice[userprobabilities -> oneuser, 88];
   (*add up how many users fell in each bin {1,2,3,4,5,6,7,8}*)
   usersamples = Total[usersperbin];
   (*Use these numbers to weigh a random choice of a single bin \
number*)
   randombin = RandomChoice[usersamples -> oneuser];
   (*For this random bin number; 
   randomly choose a corresponding number of cycles and output the \
result so it goes into the table*)   
   RandomInteger[randombin.fiveyearbins]
   (*repeat for all of the nonzero pull cases*)
   ,Length[nonzeropulllocations]];
(*assign each user one of these cycle numbers, deleting all zero \
cycle numbers*)
(*Portside*)
portsideallcycles = cyclenumbers[[1 ;; Length[portside]]];
portsidecycles = DeleteCases[portsideallcycles, 0];
(*Portrear*) 
portrearallcycles = 
  cyclenumbers[[
   Length[portside] + 1 ;; Length[portside] + Length[portrear]]];
portrearcycles = DeleteCases[portrearallcycles, 0];
(*Portbottom*)
portbottomallcycles = 
  cyclenumbers[[
   Length[portside] + Length[portrear] + 1 ;; 
    Length[portside] + Length[portrear] + Length[portbottom]]];
portbottomcycles = DeleteCases[portbottomallcycles, 0];
(*Sample the side force distribution parameters*) 
sidedist = LogNormalDistribution[\[Mu], \[Sigma]];
sideBootstrap := {\[Mu], \[Sigma]} /. 
      FindDistributionParameters[
    RandomChoice[forceside, Length[forceside]], sidedist];
sidesample = sideBootstrap;
sidemu = sidesample[[1]];
sidesigma = sidesample[[2]];
(*Sample the rear force distribution parameters*)
reardist = LogNormalDistribution[\[Mu], \[Sigma]];
rearBootstrap := {\[Mu], \[Sigma]} /. 
      FindDistributionParameters[
    RandomChoice[forcerear, Length[forcerear]], reardist];
rearsample = rearBootstrap;
rearmu = rearsample[[1]];
rearsigma = rearsample[[2]];
(*Sample the bottom force distribution parameters*)
bottomdist = LogNormalDistribution[\[Mu], \[Sigma]];
bottomBootstrap := {\[Mu], \[Sigma]} /. 
      FindDistributionParameters[
    RandomChoice[forcebottom, Length[forcebottom]], bottomdist];
bottomsample = bottomBootstrap;
bottommu = bottomsample[[1]];
bottomsigma = bottomsample[[2]];
(*Make a list of sequential forces for each nonzero force and nonzero \
cycle*)
(*sideforce*)
sideforcehistories = ParallelTable[
   RandomVariate[LogNormalDistribution[sidemu, sidesigma], 
    portsidecycles[[i]]]
   , {i, 1, Length[portsidecycles]}];
(*rearforce*)
rearforcehistories = ParallelTable[
   RandomVariate[LogNormalDistribution[rearmu, rearsigma], 
    portrearcycles[[j]]]
   , {j, 1, Length[portrearcycles]}];
(*bottomforce*)
bottomforcehistories = ParallelTable[
   RandomVariate[LogNormalDistribution[bottommu, bottomsigma], 
    portbottomcycles[[k]]], {k, 1, Length[portbottomcycles]}];
(*Combine all of the force histories into one list*)
nonzeroforcehistories = 
  Join[sideforcehistories, rearforcehistories, bottomforcehistories];
(*Restore all of the zero cycles and zero forces histories*) 
numberofzeroes = mccount - Length[nonzeroforcehistories];
allforcehistories = 
  Join[nonzeroforcehistories, ParallelTable[{0}, numberofzeroes]];
(*Find the max pull force for each nonzero force history:*)
maxstresses = Map[Max, nonzeroforcehistories];
(*The max event numbers of the nonzero force histories are the eventsamples numbers:*)
maxevents = Map[Length, nonzeroforcehistories];
(*Make a list of {maxevents,maxstresses} for all nonzero force histories:*)
maxeventsandstresses = Transpose[{maxevents, maxstresses}];
(*Count the number of {0,0} elements missing from {maxevents,maxstresses}:*)
numbermissing = mccount - Length[maxeventsandstresses];
(*make that many {0,0} elements*)
zeropairs = ParallelTable[{0, 0}, numbermissing];
(*Append these {0,0} elements to the list of maxeventsandstresses*)
fullmaxeventsandstresses = Join[maxeventsandstresses, zeropairs];
maxplot = 
  ListPlot[fullmaxeventsandstresses, PlotRange -> All, Frame -> True, 
   PlotStyle -> Red, GridLines -> Full, FrameLabel -> {"x", "y"}];
(*--------------------------------*)
(*Input data, set cycle spec, and set reliability target in %:*)
pts = fullmaxeventsandstresses;
cyclespec = 110;
rtarget = 97.50;
(*--------------------------------*)
(*Carl Woll: Presort the points by y value for speed:*)
yFunction[pts_] := 
 yFunction[Length[pts], Sequence @@ Transpose@SortBy[pts, Last]]
yFunction[len_, xpts_, ypts_][percent_, x_] := 
 Module[{xin, yin}, xin = UnitStep[x - xpts];
  If[Total[xin] < percent len, $Failed, 
   Pick[ypts, xin, 1][[Floor[len percent]]]]]
(*The yFunction constructs a yFunction object that can be used to
find the y value for a given percent and x-value:*)
yf = yFunction[pts];
(*yf[target fractional reliability, target cycle number]:*)
stressspec = yf[rtarget/100, cyclespec];
Print[{cyclespec, stressspec}];
(*Calculate the actual reliability in %*)
reliability = 
  100.*Count[
     fullmaxeventsandstresses, {u_, v_} /; 
      u < cyclespec && v < stressspec]/mccount;
Print[reliability];
(*Envelope of spec box:*)
envelope = 
  Plot[yf[reliability/100, x], {x, 0, 1400}, PlotStyle -> Blue];
(*Plot the data and the envelope:*)
Show[maxplot, envelope]
(*--------------------------------*)

Answer 1

Updated with faster method

Here's a function that takes as arguments the points, the percent, and an x-value:

y[pts_, percent_, x_] := Module[{xin, yin},
    xin = UnitStep[x - pts[[All, 1]]];
    If[Total[xin] < percent Length[pts],
        $Failed,

        yin = Pick[pts[[All, 2]], xin, 1];
        yin[[First @ Ordering[yin, {Floor[Length[pts] percent]}]]]
    ]
]

Example showing timing for a million points:

SeedRandom[1]
pts = RandomReal[1, {10^6, 2}] . {{1, 0}, {0, 5}};

y[pts, .3, .8] //AbsoluteTiming
y[pts, .7, .2] //AbsoluteTiming

{0.140796, 1.87733}

{0.027093, $Failed}

Presorting

If we presort the points by y value, we can speed this up considerably. Here is an approach that does this.

yFunction[pts_] := yFunction[Length[pts], Sequence @@ Transpose @ SortBy[pts, Last]]

yFunction[len_, xpts_, ypts_][percent_, x_] := Module[{xin, yin},
    xin = UnitStep[x - xpts];
    If[Total[xin] < percent len,
        $Failed,

        Pick[ypts, xin, 1][[Floor[len percent]]]
    ]
]

The idea is that yFunction constructs a yFunction object that can be used to find the y value for a given percent and x-value. Here is the same example:

yf = yFunction[pts];

yf[.3, .8] //AbsoluteTiming
yf[.7, .2] //AbsoluteTiming

{0.021508, 1.87733}

{0.00715, $Failed}

About 7 times faster.

Visualizations

Visualization of y vs x:

Plot[yf[.3, x], {x, 0, 1}]

Rectangle visualization for a different set of 1000 points:

SeedRandom[2]
pts = RandomReal[1, {1000, 2}] . {{5, 0}, {0, 1}};
yf = yFunction[pts];

Graphics[{
    Green,
    Point[pts],
    FaceForm[None],
    EdgeForm[Red],
    Rectangle[{0, 0}, {#, yf[.2, #]}]& /@ Range[5]
}]

Answer 2

fraction = .5
n = 1000
pts = RandomReal[{0, 1}, {n, 2}];
xmin = Sort[pts[[All, 1]]][[Ceiling[fraction n]]]
y[x_, n_] := 
 Quiet@Check[ Sort[Select[pts, #[[1]] <= x  &][[All, 2]]][[n]], 0]
Graphics[{Point@
   pts, {FaceForm[None], EdgeForm[Red], Rectangle[ {0, 0}, #]} & /@
   Select[
    Table[ {x, y[x, Ceiling[fraction n]]}, {x, 
      Subdivide[xmin, 1, 6]}], #[[2]] > 0 &]}]

Answer 3

This is a straight forward, brute force solution. It returns the lower bounds of all the bins that meet your criteria. (Assuming I've read it correctly.)

HeaneysProblem[data_, maxDataValue_, binWidth_, percentage_] :=
  Module[ {binC, R, targetN},
   binC = BinCounts[data, {0, maxDataValue, binWidth},{0, maxDataValue, binWidth}];
   R = Accumulate[Map[Accumulate,binC]];
   targetN = Floor[percentage*Length[data]];
   Position[R, _?((targetN <= # < targetN + 1) &)]*binWidth
 ]

 CheckHeaneysProblem[data_, results_] := 
    Union[Map[CheckHeaneysProblem1[data, #] &, results]]
 CheckHeaneysProblem1[data_, {rx_, ry_}] := 
    Length[Cases[data, {x_, y_} /; x <= rx && y <= ry]]/Length[data]

 RandomReal[{0, 1}, {100, 2}];
 HeaneysProblem[%, 1, 0.01, 1/10];
 ListPlot[{%, %%}]
 CheckHeaneysProblem[%%%, %%]

EDIT: Runs much faster with Henrik Schumacher's suggestion of using "accumulate" instead of summation

Answer 4

This generates some data to play with:

n = 1000000;
xa = 0.;
xb = 2.;
ya = 0.;
yb = 1.;
pts = Transpose[{RandomReal[{xa, xb}, n], RandomReal[{ya, yb}, n]}];
percentage = 10.;
threshold = Length[pts] percentage 0.01;

Now I split the whole rectangle into small reactangular bins and count the number of points in each bin with BinCounts:

xresolution = 2000;
yresolution = 1000;
counts = BinCounts[
  pts, 
  {xa, xb, (xb - xa)/xresolution}, 
  {ya, yb, (yb - ya)/yresolution}
 ];

Accumate provides an efficient way to sum up the number of pts in the rectangular region to the lower left of each bin.

accumulatedcounts = Accumulate[Accumulate /@ counts];

Now we substract the threshold and apply UnitStep and Image to obtain a black white image img. (Note that compared to the original rectangle, the image appears rotated since the first coordinate axis in images points downward). A white pixel in img stands for a bin which has more than the desired percentage of points in its lower left rectangular; a black pixel indicates less than the percentage of points to the lower left. Finally, we detect the points on the border with EdgeDetect and read off the bin indices with Position.

img = Image@UnitStep[accumulatedcounts - threshold]
edges = EdgeDetect@img
idx = Position[ImageData[edges], 1];

Finally, we merely have to extract the upper right x-y-coordinates of the bins. I do it like this:

{i, j} = Transpose[idx];
binupperrightx = Rest[Subdivide[xa, xb, xresolution]];
binupperrighty = Rest[Subdivide[xa, xb, xresolution]];
result = Transpose[{binupperrightx[[i]], binupperrightx[[j]]}];
ListPlot[result, PlotRange -> {{xa, xb}, {ya, yb}}, Frame -> True]

Because we apply rasterization, the method is not exact (we can improve upon that by increasing xresolution and yresolution), but it should suffice for every day applications.