4
$\begingroup$

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

x,y plot of nonrandom data points

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]
(*--------------------------------*)
$\endgroup$
  • $\begingroup$ For a start, you can use BinCounts to collect the number of points in each reactangle... $\endgroup$ – Henrik Schumacher Apr 17 '18 at 22:52
  • $\begingroup$ Could you please elaborate? $\endgroup$ – Michael B. Heaney Apr 18 '18 at 14:53
  • 1
    $\begingroup$ 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. $\endgroup$ – Carl Woll Apr 21 '18 at 0:08
  • $\begingroup$ @CarlWoll That works, thanks! $\endgroup$ – Michael B. Heaney Apr 23 '18 at 22:41
2
$\begingroup$

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}]

enter image description here

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]
}]

enter image description here

$\endgroup$
  • $\begingroup$ That is really helpful! $\endgroup$ – Michael B. Heaney Apr 19 '18 at 21:43
  • $\begingroup$ @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? $\endgroup$ – Carl Woll Apr 20 '18 at 2:08
  • $\begingroup$ @MichaelB.Heaney I was expecting you to evaluate those using the values of x, percent and pts that had an issue. $\endgroup$ – Carl Woll Apr 20 '18 at 23:56
  • $\begingroup$ My error: your algorithm works correctly, but I was fooled by run-to-run statistical variations. I have amended the post accordingly. Thanks again! $\endgroup$ – Michael B. Heaney Apr 21 '18 at 0:45
  • $\begingroup$ Can you please show me how to generalize your code to work on points in 3 dimensions? Thank you. $\endgroup$ – Michael B. Heaney Mar 5 at 6:06
4
$\begingroup$

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]

enter image description here

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.

$\endgroup$
5
$\begingroup$

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[%%%, %%]

enter image description here

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

$\endgroup$
  • $\begingroup$ I've tried without success to generalize this to points in 3 dimensions. Can you please show me how to do this? $\endgroup$ – Michael B. Heaney Mar 4 at 20:29
6
$\begingroup$
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 &]}]

enter image description here

$\endgroup$
  • $\begingroup$ Could you please explain what the Quiet@Check step does? $\endgroup$ – Michael B. Heaney Apr 18 '18 at 22:12
  • 1
    $\begingroup$ 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. $\endgroup$ – george2079 Apr 18 '18 at 22:16

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.