# Redistributing a set of uniformly distributed numbers to an arbitrarily defined shape

Lets say I have a random number generator that spits out uniform numbers from 0 to 1

Next, I have a shape defined by a series of vertices, like { [0, 0.4], [0.5, 0.2], [1, 0.4] }

In those vertices, the first value determines the position and second the distribution weight in that point (interpolated in between).

Now I need a magic function to redistribute the uniform random numbers into a distribution that would statistically match the shape of the vertices provided.

Essentially, if I were to run 10000 random numbers through the function, sort them into buckets like { [0 - 0.1], [0.1 - 0.2] ... } and count the buckets so that [x = x, y = count], it should produce a shape similar to if I were to use the original shape as [x, y].

To illustrate

ms paint at its finest.

Anyways, the question - I'm looking for an algorithm to build that magic function.

The shapes that I want to feed it would be based on real world statistics which often don't fit any of the neat looking distributions like gauss or normal, I've seen all kinds of wonky lines.

• You previously posted a variant of this question earlier, and got both feedback on it and received an answer: cs.stackexchange.com/q/66120/755. Don't just post a new copy of the question with some changes. Instead, you should edit the original question. – D.W. Nov 18 '16 at 16:32

## 2 Answers

Here is a simple algorithm. Suppose that the shape is given by a sequence $$(x_0,y_0), \ldots, (x_n,y_n).$$ Define $A_i = (x_i - x_{i-1}) \cdot (y_{i-1} + y_i)$ for $1 \leq i \leq n$, and let $A = \sum_{i=1}^n A_i$, $\alpha_i = A_i/A$, and $\beta_i = \sum_{j=1}^i \alpha_i$.

Sample a value $x$ uniformly on $[0,1]$, and determine $1 \leq i \leq n$ such that $\beta_{i-1} \leq x < \beta_i$. Let $c = (x-\beta_{i-1})/(\beta_i-\beta_{i-1})$.

Note that $c$ is uniform on $[0,1]$.

Here is what we have accomplished so far:

• The probability of obtaining a specific $i$ is proportional to the probability that a random point generated according to the desired density lies between $x_{i-1}$ and $x_i$.

This is because the area below the density curve between $x_{i-1}$ and $x_i$ is $(x_i - x_{i-1}) \frac{y_{i-1}+y_i}{2}$.

Consider now the function \begin{align*} \varphi_i(t) &= \int_{x_{i-1}}^t y_{i-1} + (y_i-y_{i-1}) \frac{t-x_{i-1}}{x_i-x_{i-1}} \, dt \\ &= y_{i-1} (t-x_{i-1}) + (y_i - y_{i-1})\frac{(t-x_{i-1})^2}{2(x_i-x_{i-1})} \end{align*} and its normalized cousin $$F_i(t) = \frac{\varphi_i(t)}{\varphi_i(x_i)} = \frac{2y_{i-1}}{y_{i-1}+y_i} \frac{t-x_{i-1}}{x_i-x_{i-1}} + \frac{y_i-y_{i-1}}{y_i+y_{i+1}} \left(\frac{t-x_{i-1}}{x_i-x_{i-1}}\right)^2,$$ which is the CDF of your density function restricted to the interval $[x_{i-1},x_i]$.

Output $F_i^{-1}(c)$.

Since $c \sim U(0,1)$, the resulting element has CDF $F_i$, and so has the correct distribution. It remains to explicitly invert $F_i$. Let $a = \frac{y_i-y_{i-1}}{y_i+y_{i-1}}$, $b = \frac{2y_{i-1}}{y_i+y_{i-1}}$, and $s = \frac{t-x_{i-1}}{x_i-x_{i-1}}$ (note that $0 \leq s \leq 1$). If $F(t) = c$ then $as^2 + bs = c$, and so we can calculate the output $t$ as $$s = \frac{-b + \sqrt{b^2+4ac}}{2a}, \qquad t = x_{i-1} + (x_i - x_{i-1}) s.$$ The formula for $s$ simplifies to $$s = \frac{\sqrt{(1-c) y_{i-1}^2 + c y_i^2} - y_{i-1}}{y_i-y_{i-1}}.$$ We can thus rephrase the second step as

Output $$x_{i-1} + (x_i - x_{i-1}) \frac{\sqrt{(1-c) y_{i-1}^2 + c y_i^2} - y_{i-1}}{y_i-y_{i-1}}.$$

Another option for the second step is rejection sampling. Let $A$ be the area below the rectangle bounded by $(x_{i-1},0),(x_{i-1},y_{i-1}),(x_i,y_i),(x_i,0)$, and let $R$ be the upper bounding rectangle $(x_{i-1},0),(x_{i-1},y),(x_i,y),(x_i,0)$, where $y = \max(y_{i-1},y_i)$. Generate random points $(x,y)$ in $R$ repeatedly until one of them belongs to $A$, and output its ordinate $x$. Since the area of $A$ is at least half the area of $R$, on average this requires at most two samples.

In the first step you have to determine $i$ such that $\beta_{i-1} \leq x < \beta_i$. While a naive implementation would take time $O(n)$, or $O(\log n)$ using binary search trees, you can use the alias method to reduce sampling time to $O(1)$.

You can use the inverse transform sampling method, as described in Wikipedia. See also https://stats.stackexchange.com/q/176955/2921 and https://en.wikipedia.org/wiki/Ziggurat_algorithm.

• Examples of both methods appear in my answer. – Yuval Filmus Nov 18 '16 at 21:06