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