Generate random number with specific probability distribution

Question

Well consider I have a uniform random number generator. How would I craft a function which takes as parameter this RNG, such that the distribution is following a given function?

Or more to the point: I have a a function $f(x) = 0.75-0.75x^2$ (and outside the interval $[-1, 1]$ it is $0$. This function represents the distribution I wish to gain.

Now as the integral represents the chance a number is between two boundaries, I could integrate above function between $[-1, X]$ and this value would then represent the "RNG(x)". Doing so:

$$\int_{-1}^X 0.75 -0.75x^2 dx = 2 \cdot \mathrm{RNG}(x) - 0.5$$ $\mathrm{RNG}(X)$ gives a number between $[0, 1]$ so I convert it also to the correct $$-0.25X^3 + 0.75X - 0.5 = 2 \cdot \mathrm{RNG}(x) - 0.5$$

However this is an unsolvable equation? How would I make a RNG with said distribution?

Answer 1

This is not unsolvable in the strictest sense, but a different approach may be appropriate for $f$ with support $[a,b]$ and maximum $c$:

1. Generate two independant uniform random numbers $u,v\in[0,1]$.
2. Let $x=(b-a)u+a$. If $f(x)<c v$ return $x$. Otherwise go back to step 1.

This is especially nice if $f$ differs only slightly from a uniform distribution. Alternatively, this rejection method can also be refined if one can generate $(u,v)$ uniformly in another shape that approximates the graph of $f$. The simplest refinement would be to split the ractnagle $[a,b]\times[0,c]$ into several smaller rectangles or trapezoids.

Answer 2

What about an approximation:

1. draw $N$ random numbers from intervall $[-1,1]$
2. calculate probability of each number with your $f\left(x\right)$
3. draw one on the number with probability propotional to your calculated probability.

example matlab/octave code:

allVals = zeros(1,10000);
for i = 1: 10000;
    % draw n numbers
    vals = rand(1,1000)*2-1;
    % calculate probability of each
    prob = 0.75 - 0.75.*vals.^2;
    % select proportoinal to probability
    idx = find(rand(1)*sum(prob) < cumsum(prob),1);
    allVals(i) = vals(idx);
end

figure(1);
clf();
plot(-1:0.01:1,0.75-0.75*(-1:0.01:1).^2,'b-')
figure(2);
clf();
hist(allVals,50);

Answer 3

My approach is using CDF Transformation method.

$$f(x) = 0.75 - 0.75x^2$$

$$-1\leq x \leq 1$$

add one to the range

$$0\leq x+1 \leq 2$$

Divide by 2

$$0\leq \frac{(x+1)}{2} \leq 1$$

$$P(Y\leq y) = F_{Y}(y)$$

$$P(\frac{(X+1)}{2}\leq y)$$

$$P(X\leq (2y-1))$$

$$F_{X}(x) = (.75x - .75\frac{x^3}{3})|^{x}_{-1} = (\frac{3}{4}x- \frac{1}{4} x^3) - (-.75+.75*\frac{1}{3})$$

$$F_{X}(x) = (.75x - .75\frac{x^3}{3})|^{x}_{-1} = (\frac{3}{4}x- \frac{1}{4} x^3) +\frac{1}{2}$$

$$F_{X}(2y-1) = F_{Y}(y)=(\frac{3}{4}(2y-1) - \frac{1}{4} (2y-1)^3)+\frac{1}{2}$$

where $0<y<1$

Now y is approximately a uniformily distributed variable in (0,1)

Can you take it from here.