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,
>           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,
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

```

• Prev by Date: Washington DC Area Mathematica Special Interest Group
• Next by Date: Re: Exception message from java: MathLink connection was lost.
• Previous by thread: Re: Random number with custom distribution
• Next by thread: How to lay out a grid of plots with frame labels, but no spaces