Re: random variable

• To: mathgroup at smc.vnet.net
• Subject: [mg102854] Re: [mg102833] random variable
• From: Bob Hanlon <hanlonr at cox.net>
• Date: Sat, 29 Aug 2009 06:33:49 -0400 (EDT)

```Clear[pdf, cdf, m, s, b, n1, n2];

Define some distribution such as

pdf[b_, n1_, n2_, x_] :=
(n2*b^(-n1/n2))/Gamma[n1/n2]*
x^(n1 - 1)*Exp[-x^n2/b]

Check the pdf

Assuming[{b > 0, n1 > 0, n2 > 0},
Integrate[pdf[b, n1, n2, x],
{x, 0, Infinity}]]

1

Define the cdf (this might have to be done numerically)

cdf[b_, n1_, n2_, x_] = Assuming[
{n1 > 0, n2 > 0, b > 0, x >= 0},
Simplify[
Integrate[pdf[b, n1, n2, t], {t, 0, x}]]]

1 - Gamma[n1/n2, x^n2/b]/Gamma[n1/n2]

Check the cdf

Assuming[{b > 0, n1 > 0, n2 > 0},
Limit[cdf[b, n1, n2, x], x -> #]] & /@
{0, Infinity}

{0,1}

Determine the mean

m[b_, n1_, n2_] = Assuming[
{n1 > 0, n2 > 0, b > 0},
Integrate[x*pdf[b, n1, n2, x],
{x, 0, Infinity}]]

(b^(1/n2)*Gamma[(n1 + 1)/n2])/Gamma[n1/n2]

Determine the standard deviation

s[b_, n1_, n2_] = Assuming[
{n1 > 0, n2 > 0, b > 0},
Simplify[Sqrt[
Integrate[(x - m[b, n1, n2])^2*pdf[b, n1, n2, x],
{x, 0, Infinity}]]]]

(b^(1/n2)*Sqrt[Gamma[n1/n2]*Gamma[(n1 + 2)/n2] -
Gamma[(n1 + 1)/n2]^2])/Gamma[n1/n2]

The inverse cdf is determined numerically

icdf[b_?NumericQ, n1_?NumericQ,
n2_?NumericQ, q_?NumericQ] :=
Module[{mu = m[b, n1, n2]},
Chop[x /. FindRoot[cdf[b, n1, n2, x] == q,
{x, mu}][[1]]]]

A random number drawn from the distribution is (given that it must be nonnegative for this distribution)

rand[b_?NumericQ, n1_?NumericQ,
n2_?NumericQ] := Module[{r = I},
While[Head[r] === Complex || r < 0,
r = icdf[b, n1, n2, RandomReal[]]]; r]

Defining a specific distribution

{b, n1, n2} =
Ceiling[#, .25] & /@
{2 RandomReal[], 4 RandomReal[], 3 RandomReal[]}

{0.75,4.,1.25}

Generating the data

Quiet[data = Table[rand[b, n1, n2], {400}]];

Show[{
Histogram[data, Automatic, "ProbabilityDensity"],
Plot[Tooltip[pdf[b, n1, n2, x]],
{x, 0, m[b, n1, n2] + 3 s[b, n1, n2]},
PlotStyle -> Directive[Red, Thick],
PlotRange -> All]}]

Bob Hanlon

---- omar bdair <bdairmb at yahoo.com> wrote:

=============
I want to ask, how can I generate a random vaiable from some probability density functions which are not well-known? I mean, if we have some pdf which is not normal, binomial, weibull, ... but the only thing I know that it is a log-concave function, then how can I generate a number of random variables?

```

• Prev by Date: Update on earlier post of Problems encountered with Mathematica & Snow
• Next by Date: Problem with recent file list in Mathematica 7
• Previous by thread: Re: random variable
• Next by thread: simplex representation