Re: Random number with custom distribution
- To: mathgroup at smc.vnet.net
- Subject: [mg107481] Re: Random number with custom distribution
- From: Mark Fisher <particlefilter at gmail.com>
- Date: Sun, 14 Feb 2010 05:58:31 -0500 (EST)
- References: <hl39bp$fdm$1@smc.vnet.net> <hl5uh3$ort$1@smc.vnet.net>
On Feb 13, 5:22 am, sashap <pav... at gmail.com> wrote:
> On Feb 12, 4:08 am, wiso <giurr... at gmail.com> wrote:
>
> > Dear experts, is there a simple way to generate a random sample of a
> > custom distribution? With Random I can generate random numbers
> > according to the predefined distribution (Uniform, Normal, ...) but
> > what about custom distribution (in my case a polynomial function)?
>
> On Feb 12, 4:08 am, wiso <giurr... at gmail.com> wrote:
>
> > Dear experts, is there a simple way to generate a random sample of a
> > custom distribution? With Random I can generate random numbers
> > according to the predefined distribution (Uniform, Normal, ...) but
> > what about custom distribution (in my case a polynomial function)?
>
> Since you PDF is polynomial, one approach is to
> 1. map domain of your distribution to 0...1
> 2. represent your polynomial as a linear combination
> of Sum[ w[n, m]*PDF[ BetaDistribution[ n + 1, m + 1], {n, 0,
> ord}, {m, 0, ord-n}]
> where is no less than the order of your polynomial and all w[n,m]
>= 0.
>
> 3. Then generate each random number in two steps:
> a. choose n using RandomChoice with weights w[n, m]
> b. Generate random real from BetaDistribution corresponding
> to values of {n, m} chosen in step a.
> 4. remap back to your original domain
>
> Code implementing this approach:
>
> RandomRealPolynomial[poly_, {x_, a_, b_}, len_ ] :=
> Module[{w, t, pdf, weights, coeffs, sol, k},
> pdf = ((b - a) poly /. x -> a + (b - a) t);
> k = Exponent[poly, x];
> While[True,
> coeffs = Flatten[Table[w[n, m], {n, 0, k}, {m, 0, k - n}]];
> sol = FindInstance[
> Apply[And,
> Thread[CoefficientList[
> Sum[w[n, m] (m + n + 1) Binomial[m + n, n] (1 - t)^m t^
> n, {n, 0, k}, {m, 0, k - n}] - pdf, t] == 0]] &&
> Apply[And, Thread[coeffs >= 0]], coeffs, Reals];
> If[sol =!= {},
> weights = coeffs /. First[sol];
> Break[];
> ]; k++];
> ord = RandomChoice[
> weights ->
> Flatten[Table[
> BetaDistribution[n + 1, m + 1], {n, 0, k}, {m, 0, k - n}]],
> len];
> (b - a) Table[RandomReal[dist], {dist, ord}] + a
> ]
>
> As a proof that it works try
>
> rvec = RandomRealPolynomial[3/28 (x^2 + x + 1), {x, -2, 2}, 10^6];
>
> Show[Histogram[rvec, Automatic, "ProbabilityDensity"],
> Plot[3/28 (x^2 + x + 1), {x, -2, 2}, PlotStyle -> Red]]
>
> Hope this helps,
>
> Oleksandr Pavlyk
> Wolfram Research
Oleksandr Pavlyk's solution is very nice. I hadn't realized one could
always represent a non-negative polynomial on the unit interval as a
mixture of Beta distributions.
As a small contribution, I have speeded up the random number
generation part. As originally written, the code makes multiple calls
to each of the Beta components. It's faster to make a single call to
each component and then scramble the result with RandomSample. (I've
also slightly rewritten the While statement.)
RandomRealPolynomial1[poly_, {x_, a_, b_}, len_] :=
Module[{w, t, pdf, weights, coeffs, sol, k, dists, tally},
pdf = ((b - a) poly /. x -> a + (b - a) t);
k = Exponent[poly, x] - 1;
sol = {};
While[sol === {},
k++;
coeffs = Flatten[Table[w[n, m], {n, 0, k}, {m, 0, k - n}]];
sol = FindInstance[
Apply[And,
Thread[CoefficientList[
Sum[w[n, m] (m + n + 1) Binomial[m + n, n] (1 - t)^m t^
n, {n, 0, k}, {m, 0, k - n}] - pdf, t] == 0]] &&
Apply[And, Thread[coeffs >= 0]], coeffs, Reals]
];
weights = coeffs /. First[sol];
tally = Tally[RandomChoice[weights -> Range[Length[weights]], len]];
dists =
Flatten[Table[
BetaDistribution[n + 1, m + 1], {n, 0, k}, {m, 0, k - n}]][[
tally[[All, 1]] ]];
a + (b - a) RandomSample[
Flatten[MapThread[RandomReal, {dists, tally[[All, 2]]}]]]
]
ran1 = RandomRealPolynomial1[3/28 (1 + x + x^2), {x, -2, 2},
10^6]; // Timing
On my system, this is about 5 times as fast as the original.
--Mark