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

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,
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;
(*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[];
sidesigma = sidesample[];
(*Sample the rear force distribution parameters*)
reardist = LogNormalDistribution[\[Mu], \[Sigma]];
rearBootstrap := {\[Mu], \[Sigma]} /.
FindDistributionParameters[
RandomChoice[forcerear, Length[forcerear]], reardist];
rearsample = rearBootstrap;
rearmu = rearsample[];
rearsigma = rearsample[];
(*Sample the bottom force distribution parameters*)
bottomdist = LogNormalDistribution[\[Mu], \[Sigma]];
bottomBootstrap := {\[Mu], \[Sigma]} /.
FindDistributionParameters[
RandomChoice[forcebottom, Length[forcebottom]], bottomdist];
bottomsample = bottomBootstrap;
bottommu = bottomsample[];
bottomsigma = bottomsample[];
(*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] (*--------------------------------*)  • For a start, you can use BinCounts to collect the number of points in each reactangle... Apr 17, 2018 at 22:52 • Could you please elaborate? Apr 18, 2018 at 14:53 • Use Plot[yf[reliability/100, x], {x, 0, 1400}, PlotRange->All, PlotStyle -> Blue] instead. Otherwise Plot stops plotting the function when the y value exceeds some automatically determined maximum. Apr 21, 2018 at 0:08 • @CarlWoll That works, thanks! Apr 23, 2018 at 22:41 ## 4 Answers 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
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 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 }] • That is really helpful! Apr 19, 2018 at 21:43 • @MichaelB.Heaney Do you mean it returns $Failed? An example would help. If you can't provide an example, then what does: Total[UnitStep[x - pts[[All, 1]]]] and percent Length[pts] return? Apr 20, 2018 at 2:08
• @MichaelB.Heaney I was expecting you to evaluate those using the values of x, percent and pts that had an issue. Apr 20, 2018 at 23:56
• My error: your algorithm works correctly, but I was fooled by run-to-run statistical variations. I have amended the post accordingly. Thanks again! Apr 21, 2018 at 0:45
• Can you please show me how to generalize your code to work on points in 3 dimensions? Thank you. Mar 5, 2019 at 6:06
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, #[] <= 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]}], #[] > 0 &]}] • Could you please explain what the Quiet@Check step does? Apr 18, 2018 at 22:12
• y[x,n] would throw an error if x is too small (less than xmin) It just suppresses the error and makes the function return 0. Apr 18, 2018 at 22:16

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

• I've tried without success to generalize this to points in 3 dimensions. Can you please show me how to do this? Mar 4, 2019 at 20:29

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.