How to generate a random number according to this distribution?

[wdlang]

f(x) = 1/(1+x^{3/2})

x>0

how to generate such a random number?

 

[mathman]

The most direct method needs an explicit formula for P(x) =c∫[0,x]f(u)du, where c is defined so that P(∞) = 1. The next step is find an expression for x as the inverse of P.

Now generate a random number (random variable uniform between 0 and 1) and compute x using the inverse function.

 

[wdlang]

The most direct method needs an explicit formula for P(x) =c∫[0,x]f(u)du, where c is defined so that P(∞) = 1. The next step is find an expression for x as the inverse of P.

Now generate a random number (random variable uniform between 0 and 1) and compute x using the inverse function.

but this function seems to be non-integrable

 

[digfarenough]

I might be able to help a little to get you started. Correct me if I'm wrong, but:

First, I assume you mean probability distribution, rather than a more general distribution. In that case the integral of the function must equal 1, so you need the scaling factor c that mathman mentioned.

In this case, the integral from 0 to infinity of f(x) is $$\int_0^\infty f(x) dx = 4\pi/(3\sqrt{3})$$

So your probability distribution is $$p(x) = f(x)\frac{3\sqrt{3}}{4\pi}, x\geq0.$$ (Do you mind that x=0 is allowed here?)

The integral of p(x) is a mess:
$$\int p(x) dx = \frac{\sqrt{3}}{4\pi} \left(2\sqrt{3}atan( \frac{2\sqrt{x}-1}{\sqrt{3}}) - 2log(1+\sqrt{x}) + log(1-\sqrt{x}+x) \right) + C$$
Use your favorite calculus software to get this instead of trying by hand! E.g. wolfram alpha

So to get the cdf you need to find $$P(x) = \int_0^x p(u) du = \frac{\sqrt{3}}{4\pi} \left(2\sqrt{3}atan( \frac{2\sqrt{x}-1}{\sqrt{3}}) - 2log(1+\sqrt{x}) + log(1-\sqrt{x}+x) \right) + \frac{1}{4}$$

...if my quick math is right. You'll want to recheck that (recall sqrt(0) = 0, log(1) = 0, atan(-1/sqrt(3)) = -pi/6, and then it's just a bunch of cancellations).

Now the hard part: You want to invert P(x), some $$g(x) = P^{-1}(x).$$ I don't know if there's an easy way of doing it. You could use a look-up table by just evaluating P(x) over a set of points and inverting that.

Then, again as mathman said, you generate numbers from U[0,1] and use g(U[0,1]) to sample from p(x). If you do this, make sure you generate a big batch of them and compare them to the original distribution to make sure no bugs or typos have been introduced.

I may have made one or more mistakes there, but this is the general way of sampling from an arbitrary distribution: calculate the cdf of your arbitrary distribution, then invert that cdf and use it to translate from U[0,1] to your desired distribution.

Hope that is slightly helpful!

 

[haruspex]

An alternative approach is to keep generating X values and challenging them: you take another random value from a uniform distribution and compare it against a threshold that depends on f(x).
That's not too hard when the range of X is finite; you generate the values uniformly across the range. With unlimited X it's trickier, but it can be done.

 

[wdlang]

I might be able to help a little to get you started. Correct me if I'm wrong, but:

First, I assume you mean probability distribution, rather than a more general distribution. In that case the integral of the function must equal 1, so you need the scaling factor c that mathman mentioned.

In this case, the integral from 0 to infinity of f(x) is $$\int_0^\infty f(x) dx = 4\pi/(3\sqrt{3})$$

So your probability distribution is $$p(x) = f(x)\frac{3\sqrt{3}}{4\pi}, x\geq0.$$ (Do you mind that x=0 is allowed here?)

The integral of p(x) is a mess:
$$\int p(x) dx = \frac{\sqrt{3}}{4\pi} \left(2\sqrt{3}atan( \frac{2\sqrt{x}-1}{\sqrt{3}}) - 2log(1+\sqrt{x}) + log(1-\sqrt{x}+x) \right) + C$$
Use your favorite calculus software to get this instead of trying by hand! E.g. wolfram alpha

So to get the cdf you need to find $$P(x) = \int_0^x p(u) du = \frac{\sqrt{3}}{4\pi} \left(2\sqrt{3}atan( \frac{2\sqrt{x}-1}{\sqrt{3}}) - 2log(1+\sqrt{x}) + log(1-\sqrt{x}+x) \right) + \frac{1}{4}$$

...if my quick math is right. You'll want to recheck that (recall sqrt(0) = 0, log(1) = 0, atan(-1/sqrt(3)) = -pi/6, and then it's just a bunch of cancellations).

Now the hard part: You want to invert P(x), some $$g(x) = P^{-1}(x).$$ I don't know if there's an easy way of doing it. You could use a look-up table by just evaluating P(x) over a set of points and inverting that.

Then, again as mathman said, you generate numbers from U[0,1] and use g(U[0,1]) to sample from p(x). If you do this, make sure you generate a big batch of them and compare them to the original distribution to make sure no bugs or typos have been introduced.

I may have made one or more mistakes there, but this is the general way of sampling from an arbitrary distribution: calculate the cdf of your arbitrary distribution, then invert that cdf and use it to translate from U[0,1] to your desired distribution.

Hope that is slightly helpful!

this is very helpful.

the inversion is not a big problem.