How can I randomly select 1000 points uniformly from the shaded area below the plotted curve?
Plot[1/π Cos[θ]^2, {θ, 0, 2 π}, Filling -> Bottom] 
$\begingroup$ $\endgroup$
2 
- 4$\begingroup$ One simple way is to sample points uniformly both above and below the curve, and keep only the latter. $\endgroup$bbgodfrey– bbgodfrey2016-11-22 04:33:28 +00:00Commented Nov 22, 2016 at 4:33
- 2$\begingroup$ Given the symmetries, if (x,y) is above the curve you can replace it with ((x+π/2) mod 2π, 1/π - y). $\endgroup$Anton Sherwood– Anton Sherwood2016-11-24 08:27:50 +00:00Commented Nov 24, 2016 at 8:27
Add a comment |
4 Answers
$\begingroup$ $\endgroup$
As noted in my comment, one approach is as follows. First, generate thousands of pairs of random numbers in the range {0, 2 π}, {0, 1/π}. Then select the first 1000 that lie below the curve.
lst = Transpose@{RandomReal[{0, 2 π}, 4000], RandomReal[{0, 1/π}, 4000]}; listsel = Select[lst, #[[2]] < 1/π Cos[#[[1]]]^2 &, 1000]; Show[Plot[1/π Cos[θ]^2, {θ, 0, 2 π}, Filling -> Bottom], ListPlot[listsel]] 
This simple process works well provided the portion of points selected is a reasonable fraction of the total number of points, as it is here.
$\begingroup$
$\endgroup$
4 There is no need in filtering out random points in a rectangle that don't fall in the prescribed region. The sampling within a region can be done directly with RandomPoint.
Specify the region:
reg = ImplicitRegion[0 <= x <= 2 Pi && 0 <= y <= 1/Pi Cos[x]^2, {x, y}] and then one can sample a point with
RandomPoint[reg] e.g., {0.39486, 0.0422331}
or several points n with RandomPoint[reg, n]. There's a warning about an unbounded region, so to keep it clean one can add bounds as a third argument to RegionPlot:
data = RandomPoint[reg, 1000, {{0, 2 Pi}, {0, 1/Pi}}]; Show[RegionPlot[reg], ListPlot[data, Frame -> True], AspectRatio -> 1/GoldenRatio] 
EDIT as per Trilarion's comment:
How does RandomPoint work internally is beyond my knowledge, but the timing analysis shows that it does not sample a rectangle and throw away the points that don't fall in the region (and even if it does, it's a way faster implementation):
reg = ImplicitRegion[0 <= x <= 10 && y >= x && y <= x + 1, {x, y}] 
n = 110000; (lst = Transpose@{RandomReal[{0, 10}, n], RandomReal[{0, 11}, n]}; sel = Select[lst, #[[1]] <= #[[2]] <= #[[1]] + 1 &, UpTo[10000]];) // AbsoluteTiming Length@sel {0.221189, Null}
9916
10000 points weren't even generated.
RandomPoint approach:
pts = RandomPoint[reg, 10000, {{0, 10}, {0, 11}}]; // AbsoluteTiming {0.049927, Null}
More than 4x faster, and all 10000 desired points are obviously generated.
- $\begingroup$ "There is no need in filtering out random points in a rectangle that don't fall in the prescribed region. The sampling within a region can be done directly with RandomPoint." Isn't it kind of doing the same as filtering out random points internally? Or how does RandomPoint work? $\endgroup$NoDataDumpNoContribution– NoDataDumpNoContribution2016-11-23 07:58:01 +00:00Commented Nov 23, 2016 at 7:58
- $\begingroup$ perhaps, it first parametrizes the implicit region in a [0,1]x[0,1] topologically equivalent object and then does non-uniform sampling. This should be faster, but you need to take into account for the differential area change between the original region and the parametrization. $\endgroup$Fabio– Fabio2016-11-23 11:31:05 +00:00Commented Nov 23, 2016 at 11:31
- 1$\begingroup$ @Trilarion, one does not necessarily have to do rejection sampling; see this for instance. $\endgroup$J. M.'s missing motivation– J. M.'s missing motivation2017-01-06 18:26:43 +00:00Commented Jan 6, 2017 at 18:26
- $\begingroup$ @J.M. Thanks for the link. The answers there are great. $\endgroup$NoDataDumpNoContribution– NoDataDumpNoContribution2017-01-07 10:28:23 +00:00Commented Jan 7, 2017 at 10:28
$\begingroup$ $\endgroup$
More of a first principles approach, use the function as a PDF to generate random x data, then for each x choose a uniformly distributed point on the vertical line {x, f[x]}:
f[x_] := 1/π Cos[x]^2 z = Integrate[f[x], {x, 0, 2 π}]; (*can use NIntegrate here if needed*) Plot[f[x], {x, 0, 2 π}, Epilog -> Point[ {#, First@RandomVariate[UniformDistribution[{0, f[#]}], 1]} & /@ RandomVariate[ProbabilityDistribution[f[x]/z, {x, 0, 2 π}], 1000] ] ] 
This is likely the fastest approach.
I think newer versions of ProbabilityDistribution may do that normalization (/z) automatically, btw.
$\begingroup$ $\endgroup$
2 In case where only the plot is given:
plot = Plot[1/π Cos[θ]^2, {θ, 0, 2 π}, Filling -> Bottom] polygons = Cases[plot // Normal, _Polygon, ∞] region = RegionUnion @@ polygons; pts = RandomPoint[region, 100]; (*quite slow*) Show[plot, Graphics@Point@pts] 
- $\begingroup$ Interesting application of
Normal. PerhapsNormalis the answer the question regarding an alternative for FullGraphics ( mathematica.stackexchange.com/questions/83648/…). $\endgroup$LouisB– LouisB2016-11-23 00:41:45 +00:00Commented Nov 23, 2016 at 0:41 - 1$\begingroup$ @LouisB Probably not, it only converts
GraphicsComplexto regular primitives. $\endgroup$Kuba– Kuba2016-11-23 06:41:48 +00:00Commented Nov 23, 2016 at 6:41