3
$\begingroup$

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.453*Exp[-1.036*x]*Sinh[Sqrt[2.29*x]]. 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.453*Exp[-1.036*x]*Sinh[Sqrt[2.29*x]],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.

Improve this question
$\endgroup$
3
  • 3
    $\begingroup$ Don't program mathematica like python... Vectorize things, for example, generating randoms one at a time is grossly inefficient. $\endgroup$
    – ciao
    Commented Nov 14, 2020 at 11:06
  • 2
    $\begingroup$ see the answers in RandomVariate from 2-dimensional probability distribution and Acceptance Rejection Method $\endgroup$
    – kglr
    Commented Nov 14, 2020 at 11:33
  • 2
    $\begingroup$ 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. $\endgroup$
    – b3m2a1
    Commented Nov 14, 2020 at 20:03

1 Answer 1

Reset to default
5
$\begingroup$

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
*)
Improve this answer
$\endgroup$
9
  • $\begingroup$ 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. $\endgroup$
    – Michael E2
    Commented Nov 14, 2020 at 20:39
  • $\begingroup$ It is nice code (+1). It is not so clear for schoolboy how do you calculate sumAR2? $\endgroup$
    – Alex Trounev
    Commented Nov 16, 2020 at 20:11
  • $\begingroup$ 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! $\endgroup$
    – jojell
    Commented Nov 17, 2020 at 17:13
  • $\begingroup$ @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. $\endgroup$
    – Michael E2
    Commented Nov 17, 2020 at 18:31
  • $\begingroup$ @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. $\endgroup$
    – Michael E2
    Commented Nov 17, 2020 at 18:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.