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