Mathematica Asked on July 29, 2021
I need to generate a random number (real) and then test if this number works on my pdf. if it works i append it to a list. if it doesn’t, i reject it.
Until now i have done this:
f[x_] := x/2; "with 0<=x<=2"
this is the function i’m using
testexc = {};
xc = {};
These are the lists where i want to append numbers
While[Length[testexc] <= 10000,
AppendTo[testexc, RandomReal[{0, 1}]]];
This is how i generate 10000 Random Numbers within the specified range
What i want now is: I need 10000 random numbers that works with my function. And i couldnt find out any way of doing it.
Here, try this; it should be faster:
pdf = x [Function] x/2;
cdf = x [Function] Evaluate[Integrate[pdf[t], {t, 0, x}]]
cdfinv = y [Function] 2 Sqrt[y]
rand = cdfinv[RandomReal[{0, 1}, {1000000}]]; // RepeatedTiming // First
0.009
One million random number is a percent of a second.
Plotting a histogram to check the distribution is correct:
Histogram[rand]
Actually, there is also a built-in method for this. It goes like this:
distro = ProbabilityDistribution[x/2, {x, 0, 2}];
rand = RandomVariate[distro, 1000000]; // RepeatedTiming // First
0.08
For some reason, it is significantly slower...
If you insist on using the "acceptance/rejection method" (better know as Monte Carlo method, you can do this:
n = 2000000;
First@RepeatedTiming[
x = RandomReal[{0, 2}, n];
y = RandomReal[{0, 1}, n];
rand = Pick[x, UnitStep[Subtract[pdf[x], y]], 1];
]
0.062
This generates about a million random numbers with 0.062 seconds. I would strongly discourage methods that use Append
repeatedly, because they will have quadratic complexity and be very memory bound (each time Append
is called, you have to copy the full array).
Internal`Bag
This is very, very slow, also because random numbers a more efficiently created in bulks instead of one-by-one.
n = 1000000;
Do[
x = RandomReal[{0, 2}];
y = RandomReal[{0, 1}];
If[y <= pdf[x], Internal`StuffBag[bag, x]];
,
{n}
]; // RepeatedTiming // First
rand = Internal`BagPart[bag, All];
3.2
This takes about 3.2 seconds...
Compile
and Internal`Bag
Compiling the latter can be faster by more than two orders of magnitude, though.
cf = Block[{x},
With[{pdfx = pdf[x]},
Compile[{{n, _Integer}},
Block[{x, y, bag},
bag = Internal`Bag[Most[{0.}]];
Do[
x = RandomReal[{0, 2}];
y = RandomReal[{0, 1}];
If[y <= pdfx, Internal`StuffBag[bag, x]];
,
{n}
];
Internal`BagPart[bag, All]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
]
]
];
n = 1000000;
rand = Join @@ cf[ConstantArray[n/4, {4}]]; // RepeatedTiming // First
0.022
Answered by Henrik Schumacher on July 29, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP