       Re: Random number with custom distribution

• To: mathgroup at smc.vnet.net
• Subject: [mg107457] Re: Random number with custom distribution
• From: sashap <pavlyk at gmail.com>
• Date: Sat, 13 Feb 2010 05:22:35 -0500 (EST)
• References: <hl39bp\$fdm\$1@smc.vnet.net>

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

```

• Prev by Date: Re: Re: Normality test
• Next by Date: Re: How to scope PopupMenu values in a DynamicModule?
• Previous by thread: Re: Random number with custom distribution
• Next by thread: Re: Random number with custom distribution