Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

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)
  • Reply-to: hanlonr at cox.net

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