# Distributing values according to a given continuous probability distribution

I have a series of continuous probability distributions and the intervals for each distribution of possible values. Now I want to use the distributions to distribute values of a random variable among agents. I have thought of one way to do this, and I'm curious if this is an acceptable method. This is how I would do it:

Let P(A<X<B) be a given pdf, where A,B are min and max values for random
variable X.
Divide the spread between A and B into an arbitrarily large amount of 'chunks'
or ranges.
For each chunk, represented as a range, say between a and b, make a table
where P(a<X<b) relates to the range a,b.
For each agent receiving a value of random variable X, use a random number
generator, look up the result in the table, give the agent a value in the
range (a,b).


Assuming that this is an acceptable technique, I'm still curious as to what an acceptable amount of granularity for the 'chunks' would be. I hope this wan't too confusing and any feedback would be much appreciated.

edit: All of my functions take the form of f(x) = a/(b*x)+c over some interval. They are all easily integrable.

-
It rather depends on the smoothness of your pdf as you are going to interpolate. But if you can calculate $\Pr (a<X<b)$ then presumably you have some way of integrating the pdf, which rather changes possible approaches to the question. Perhaps you could tell us what form the pdf takes. –  Henry Jul 21 '11 at 20:28
If you can easily integrate the pdf, to acquire the cdf, can you easily invert the CDF, $F^{-1}$? If so, just plug uniform random numbers into $F^{-1}$ and the resulting numbers will be a draw from your distribution –  Macro Jul 21 '11 at 21:26
@Macro can you elaborate on this, please? If I understand you correctly, this would be using the cdf where P(X<B) for some value B, the lower bound always being A. Then I would solve for B, and have the invers function? –  Graham Jul 21 '11 at 21:40
Let $F(x) = P(X \leq x)$. Then set $y = F(x)$ and solve for $x$. The general solution to that inversion problem is the inverse CDF, $F^{-1}$. If you generate Uniform(0,1) random numbers and plug them into $F^{-1}$ you have random numbers with CDF $F$. For example, $F(x) = 1 - e^{-x}$ for $x > 0$ is the exponential CDF. The general solution to $y = 1 - e^{-x}$ in terms of $x$ is $x = F^{-1}(y) = -\log(1-y)$. Then, it follows that if $u$ is a uniform(0,1) random variable, then $-\log(1-u)$ has an exponential distribution. –  Macro Jul 21 '11 at 21:46
@Macro - I think that is an answer rather than a comment, and was what I was hinting at. My only slight extension is that even if you do not have a closed form for $F^{-1}$, you may be able to find a reasonable approximation. –  Henry Jul 22 '11 at 6:58
show 1 more comment

In your problem it sounds like the CDF will have a simple closed form that is easily inverted. In that case I'd suggest what is called the "inversion method" for generating from the distribution:

Let $F(x)=P(X≤x)$. Then set $y=F(x)$ and solve for $x$. The general solution to that inversion problem is the inverse CDF, $F^{-1}$. If you generate Uniform(0,1) random numbers and plug them into $F^{-1}$ you have random numbers with CDF $F$.

For example, $F(x)=1−e^{-x}$ for $x>0$ is the exponential CDF. The general solution to $y=1−e^{−x}$ in terms of $x$ is

$$x=F^{-1}(y)=−\log(1−y).$$ Then, it follows that if $u$ is a ${\rm Uniform}(0,1)$ random variable, then $-\log(1−u)$ has an exponential distribution.

Edit: From your problem statement, the density is proportional to $f(x) = \frac{a}{bx} + c$ on some interval $(A,B)$. Therefore, the CDF is

$$F(x) = \frac{1}{k} \int_{A}^{x} \frac{a}{by} + c \ dy = \frac{1}{k} \cdot \left( \frac{a}{b} \cdot \log(x/A) + c(x - A) \right)$$

where $k$ is the normalizing constant, which is calculated here. This closed form of the inverse of this CDF involves the Lambert W Function (thanks, @whuber), so I'm going to be lazy and solve for it numerically. Here is a very crude method of doing that in R using bisection, implemented as the uniroot function:

# Find the point at which the CDF above equals p, for parameters values
# a,b,c and predefined domain (A,B)
F_inv = function(p,a,b,c,A,B)
{
k = (a * log(B/A))/b + c(B-A)
br = c(A,B)
G = function(x) (a * log(x/A)/b + c(x-A))/k - p
return( uniroot(G,br)\$root )
}

# x is a sample from this distribution with a=2, b=5, c=.8, A=10, B=20.
x = rep(0,1000)
for(i in 1:1000) x[i] = F_inv(runif(1),2,5,.8,10,20)

-
+1 This CDF has a "closed form inverse" using the "ProductLog" (Lambert W function). Wolfram Alpha will compute it and lists some of its properties. –  whuber Jun 4 '12 at 14:34
Thanks @whuber, I didn't know that. I wrote it down and stared at it for a minute and it looked unsolvable. But, it seems like there's a special function for many such problems. –  Macro Jun 4 '12 at 14:43
There are R implementations of W: here's one promising package (I haven't checked it out) at cran.r-project.org/web/packages/LambertW/index.html and another is blogged at r-bloggers.com/…. –  whuber Jun 4 '12 at 14:48