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?