How can I make my accept-reject code run faster?

Question

For a school project we were asked to do a accept-reject method, which can calculate the mean for functions that are not very easy to solve analytically. This I did in Mathematica 12.1 and it worked, however quite slow. If i code the same principle in python it takes 6 times less time, making the ego of my python using friend way too big. Is there a way to make the following "Acceptance-Rejection" code run significantly faster?

(*Initiating a random seed for reproducibilty together with some parameters*)
ceiling = 0.3582061533158414; w=20.; h = 1.05*ceiling;
SeedRandom[4242424242];
(the function ar will generate 2 random points x and y in the perimeter of a triangle enclosing 
f[x] = 0.453Exp[-1.036*x]Sinh[Sqrt[2.29x]]. if y<f[x], the point falls in the graph and x is
returned. if not, the function will run itself again until a x,y pair is found that is
enclosed by f[x])
ar :=  (
y = h(1.-Sqrt[RandomReal[]]);
x = RandomReal[{0.,w*(h-y)/h}];
If[y<=0.453Exp[-1.036x]Sinh[Sqrt[2.29x]],Return[x]];
ar
)
(generate n points using the accept reject method and add this value ei to the sumAR. also add ei^2 to the sumAR2 for calculating the variance.)
n = 10^6; sumAR = 0; sumAR2 = 0; x = 0; y = 0;
Timing[
Do[
ei = ar; sumAR += ei; sumAR2 += ei^2
,n]
]
mean = sumAR/n
variance = sumAR2/n - mean^2
deviation = Sqrt[variance]

PS: I first had f[x] defined and then would run y<=f[x] since it is easier to read, but not having to call f[x] saves around 2 seconds in total and i only use f[x] once so yea.

Don't program mathematica like python... Vectorize things, for example, generating randoms one at a time is grossly inefficient. — ciao, Nov 14 '20 at 11:06
see the answers in RandomVariate from 2-dimensional probability distribution and Acceptance Rejection Method — kglr, Nov 14 '20 at 11:33
Even in python, this kind of approach will be super slow. All modern numerical software packages (Mathematica/Julia/Python+NumPy) work much faster when allowed to use vectorized operations. — b3m2a1, Nov 14 '20 at 20:03

Michael E2 · Answer 1 · 2020-11-17T20:34:51.973

5

An example of what the comments are talking about (OP's code takes 25-26 sec. versus 0.32 sec. below):

nn = 10^6;
pp = 0;         (* points generated so far *)
sumAR = 0; sumAR2 = 0;
aratio = 1.;    (* area proportion estimate *)
While[pp < nn,
  yy = h*(1. - Sqrt[RandomReal[1, (nn - pp)/aratio // Ceiling]]);
  xx = RandomReal[1, (nn - pp)/aratio // Ceiling]*w*(h - yy)/h;
  inOut = UnitStep[      (* inside=1, outside=0 *)
    0.453*Exp[-1.036*xx]*Sinh[Sqrt[2.29*xx]] - yy
    ];
  {sumAR, sumAR2} += Total /@ {#, #^2} &@Take[
    Pick[xx, inOut, 1],  (* accept=1, reject=0 *)
    UpTo[nn - pp]];
  pp += If[pp == 0,      (* set area ratio on first step *)
      aratio = N@#/nn; #, #] &@Total@inOut;
  ] // AbsoluteTiming
mean = sumAR/nn
variance = sumAR2/nn - mean^2
deviation = Sqrt[variance]

(*
{0.320782, Null}
1.98308
2.42222
1.55635
*)

edited Nov 17 '20 at 20:34

answered Nov 14 '20 at 20:11

Michael E2

235,386
17
334
747

For simplicity, there's region = ImplicitRegion[0 <= x <= 20 && 0 <= y <= 0.453*Exp[-1.036*x]*Sinh[Sqrt[2.29*x]] && 0 <= y <= h, {x, y}]; RandomPoint[DiscretizeRegion[region, MaxCellMeasure -> {"Length" -> 0.1}], 10^6][[All, 1]] // Mean // RepeatedTiming. DiscretizeRegion introduces an additional truncation error, and, of course, it's not the accept/reject method. It's also nearly twice as slow as the one in my answer. – Michael E2 Nov 14 '20 at 20:39
It is nice code (+1). It is not so clear for schoolboy how do you calculate sumAR2? – Alex Trounev Nov 16 '20 at 20:11
Thank you for the reply and great upgrade of my code! I had to implement the sumAR2 myself but the structure of the code is very interesting. A bit hard to understand completely for a newbie like me. Also to me it is somewhat counte- intuitive that working with a list over a million values long is quicker than just use 1 variable to sum over, but case in point i guess. Again, thank you! – jojell Nov 17 '20 at 17:13
@jojell You're welcome! The intuition comes from knowing roughly how current CPUs work (memory management, pipelining, FP units in the CPU) that allow an FP arithmetic operation that individually takes say 100 clock cycles to be done on average of more than one per cycle when the operation is performed on a long list of floating-point numbers. – Michael E2 Nov 17 '20 at 18:31
@AlexTrounev I thought the premise of Wolfram is that computer-based math is perfect for a schoolboy (well, I'm sure they would say "schoolchildren"). I assume it was a rhetorical question, but to get sumAR2, starting with the x coordinates of the points we Pick-ed that are inside the graph and Take-ing up to the desired number nn of points, we add up the Total of the x coordinates. – Michael E2 Nov 17 '20 at 18:40
@MichaelE2 Can you add to your code something like sumAR2=0; and sumAR2 += Total@Take[Pick[xx^2, inOut, 1], UpTo[nn - pp]]; for the future generations? :) – Alex Trounev Nov 17 '20 at 20:18
@AlexTrounev Oops, I totally missed that there were two sum variables. Thanks for persisting! – Michael E2 Nov 17 '20 at 20:39
@jojell Sorry, I lost sight of there being a second sumAR2 variable. Fixed now. – Michael E2 Nov 17 '20 at 20:41
@MichaelE2 Ok! Now it looks perfect. – Alex Trounev Nov 17 '20 at 20:48

How can I make my accept-reject code run faster?

1 Answers1