TransWikia.com

Accelerate the speed of selection in Random

Mathematica Asked by kile on August 4, 2020

Tossing 4 points on a circle at a time, and calculate the chance of all 4 points on the upper part of the circle.

An illustration for this

This is very simple by math:

$$left( frac{1}{2} right)^4 = 0.0625$$

But I want do a simulation in Mathematica.

po = RandomVariate[UniformDistribution[{{0, 1}, {0, 2 Pi}}], 40000000];
f[m_, n_] := {m Cos[n], m Sin[n]};
Apply[f, po, {1}] // Partition[#, 4] & // 
   Select[#, AllTrue[Last /* GreaterEqualThan[0]]] & // 
  Length // AbsoluteTiming

it takes around 140 seconds for this code to finish.

Is there a way to speed this up?

2 Answers

We get additional speed-up combining RandomPoint with UnitStep,Total and Min or Count:

n = 10^6;

SeedRandom[1];

RepeatedTiming[
 Total[Min /@ UnitStep[RandomPoint[Disk[], {n/4, 4}][[All, All, 2]]]]]
{0.099, 15733}
SeedRandom[1];

RepeatedTiming[
  Count[4] @ Total[UnitStep[RandomPoint[Disk[], {n/4, 4}][[All, All, 2]]], {2}]]
{0.11, 15733}

versus the method from Bob Hanlon's answer:

SeedRandom[1];

RepeatedTiming[
 RandomPoint[Disk[], n] // Partition[#, 4] & // 
   Select[#, AllTrue[Last /* GreaterEqualThan[0]]] & // Length]
{0.8295, 15733}

Correct answer by kglr on August 4, 2020

It is more efficient to use RandomPoint and RandomPoint distributes the points uniformly over the region.

Clear["Global`*"]

f[m_, n_] := {m Cos[n], m Sin[n]};

n = 80000;

To accurately compare the two methods, the point generation needs to be included in the timing.

RepeatedTiming[
 po = RandomVariate[
   UniformDistribution[{{0, 1}, {0, 2 Pi}}], n];
 (Apply[f, po, {1}] // Partition[#, 4] & // 
      Select[#, AllTrue[Last /* GreaterEqualThan[0]]] & // Length)/(n/4) // N]

(* {0.160, 0.06455} *)

Note that the points cluster near the origin.

ListPlot[f @@@ po, AspectRatio -> 1]

enter image description here

RepeatedTiming[
 ((pp = RandomPoint[Disk[], n]) // Partition[#, 4] & // 
      Select[#, AllTrue[Last /* GreaterEqualThan[0]]] & // Length)/(n/4) // N]

(* {0.0423, 0.0619} *)

ListPlot[pp, AspectRatio -> 1]

enter image description here

Answered by Bob Hanlon on August 4, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP