Generating Random Samples from Non-Identical Truncated Order Statistics

TL;DR

This paper introduces an efficient algorithm for sampling from the bounded kth order statistic of independent, non-identically distributed variables, especially effective in sparse density regions with boundary constraints.

Contribution

The paper presents a novel algorithm for generating samples from non-identical truncated order statistics, improving efficiency over rejection sampling in sparse density scenarios.

Findings

Algorithm is faster in sparse density regions.

Effective for tight boundary conditions.

Outperforms rejection sampling in specific cases.

Abstract

We provide an efficient algorithm to generate random samples from the bounded kth order statistic in a sample of independent, but not necessarily identically distributed, random variables. The bounds can be upper or lower bounds and need only hold on the kth order statistic. Furthermore, we require access to the inverse CDF for each statistic in the ordered sample. The algorithm is slightly slower than rejection sampling when the density of the bounded statistic is large, however, it is significantly faster when the bounded density becomes sparse. We provide a practical example and a simulation that shows the superiority of this method for sparse regions arising from tight boundary conditions and/or over regions of low probability density.

Equations

$$C^{\prime }=([a_1,1]\times [a_2,1])\setminus ([b_1,1]\times [b_2,1])$$

$$R_{s}^{\prime}=\left\{(u_{1}^{\prime},u_{2}^{\prime}):\begin{array}[]{cc}\forall j\in s,&a_{j}<u_{j}^{\prime}<b_{j}\\ \forall j\notin s,&b_{j}<u_{j}^{\prime}<1\end{array}\right\}$$

$$R_s=\{(u_1,u_2):0<u_1<1;\sum _{j=1}^{s-1}{f}_j<u_2<\sum _{j=1}^s{f}_j\}$$

$$\displaystyle g_{s}(u_{j})=\left\{\begin{array}[]{cc}(b_{j}-a_{j})u_{j}+a_{j}&\text{if }j\in s,j=1\\ \frac{b_{j}-a_{j}}{f_{s}}(u_{j}-\sum\limits_{l=1}^{s-1}f_{l})+a_{j}&\text{if }j\in s,j=2\\ (1-b_{j})u_{j}+b_{j}&\text{if }j\notin s,j=1\\ \frac{1-b_{j}}{f_{s}}(u_{j}-\sum\limits_{l=1}^{s-1}f_{l})+b_{j}&\text{if }j\notin s,j=2\end{array}\right.$$

$$C^{\prime}=\left\{(u_{1}^{\prime},\cdots,u_{N}^{\prime}):\begin{array}[]{cc}a_{j}<u_{j}^{\prime}<b_{j}&\text{for at least }k\text{ of the }j\in\{1,\cdots,N\}\\ 0<u_{j}^{\prime}<1&\text{otherwise}\end{array}\right\}$$

$$S_A$$

$$S_B$$

$$S=\{s=(s_A,s_B):s_A\in S_A,s_B\in S_B,s_A\cup s_B=I,s_A\cap s_B\ne \varnothing \}$$

$$R_{s=(s_{A},s_{B})}^{\prime}=\left\{(u_{1}^{\prime},\cdots,u_{N}^{\prime}):\begin{array}[]{cc}a_{j}<u_{j}^{\prime}<b_{j}&\text{if }j\in s_{A}\cap s_{B}\\ b_{j}<u_{j}^{\prime}<1&\text{if }j\in s_{A}\setminus s_{B}\\ 0<u_{j}^{\prime}<a_{j}&\text{if }j\in s_{B}\setminus s_{A}\end{array}\right\}.$$

$$V_{R_s^{\prime }}=\prod _{l\in s_A\cap s_B}(b_l-a_l)\prod _{m\in s_A\setminus s_B}(1-b_m)\prod _{n\in s_B\setminus s_A}{a}_n$$

$$R_{s}=\left\{(u_{1},\cdots,u_{N}):\begin{array}[]{cc}0<u_{j}<1&\text{if }j<N\\ \sum\limits_{l=1}^{s-1}f_{l}<u_{j}<\sum\limits_{l=1}^{s}f_{l}&\text{if }j=N\end{array}\right\}$$

$$\displaystyle u_{j}^{\prime}=g_{s}(u_{j})=\left\{\begin{array}[]{cc}(b_{j}-a_{j})u_{j}+a_{j}&\text{if }j\in s_{A}\cap s_{B},j<N\\ \frac{b_{j}-a_{j}}{f_{s}}(u_{j}-\sum\limits_{l=1}^{s-1}f_{l})+a_{j}&\text{if }j\in s_{A}\cap s_{B},j=N\\ (1-b_{j})u_{j}+b_{j}&\text{if }j\in s_{A}\setminus s_{B},j<N\\ \frac{1-b_{j}}{f_{s}}(u_{j}-\sum\limits_{l=1}^{s-1}f_{l})+b_{j}&\text{if }j\in s_{A}\setminus s_{B},j=N\\ a_{j}u_{j}&\text{if }j\in s_{B}\setminus s_{A},j<N\\ \frac{a_{j}}{f_{s}}(u_{j}-\sum\limits_{l=1}^{s-1}f_{l})&\text{if }j\in s_{B}\setminus s_{A},j=N\end{array}\right.$$

Full text

1 Introduction

We consider the problem of how to efficiently draw pseudo-random variates from a truncated distribution where the distribution is equivalent to the $k$-th order statistic of a set of $N$ independent random variables. The CDF and inverse CDF for each of the $N$ distributions is known. The main contribution of this paper is an algorithm that bounds the $N$-dimensional hypercube so that the inverse CDF sampling method may be used (see Devroye [2006]).

Consider $N$ independent random variables $X_{i}$ whose CDFs, $F_{X_{i}}$, and inverse CDFs, $F_{X_{i}}^{-1}$, are known, and we want to draw a random sample from their $k$-th order statistic, $Y=X_{(k)}$ satisfying bounds of $A<Y<B$. As the bounds narrow, rejection sampling struggles to find acceptable samples. Instead of waiting for acceptable samples, we map the samples from the $N$-dimensional unit hypercube onto a restricted subset of another $N$-dimensional unit hypercube that is normalized to represent the restricted probability space. We then apply the inverse CDF to obtain our draws from $X_{i}$.

The procedure is useful in a range of applications that support the factorization of a distribution function into $N$ known distribution functions not necessarily identical. A frequent example is drawing bounded variates from the minimum or maximum of a sample of order statistics, which arise when modeling censored or truncated data. The 2-dimensional minimum or maximum case is presented as a practical motivating example prior to the general bounded $k$-th order statistic.

2 Two-Dimensional Case

Define $Y=\min(X_{1},X_{2})$ where each $X_{i}$ are independent, and $Y$ is bounded below by $A$ and above by $B$. Further, assume we have access to their CDFs, $F_{X_{i}}$, and inverse CDFs, $F_{X_{i}}^{-1}$. Without the bound restriction, the sample is easily generated by drawing two independent uniform draws from $[0,1]$, applying the inverse CDF to each draw and returning the minimum of the two. To satisfy the boundary condition, we restrict the region as

[TABLE]

where $a_{i}=F_{X_{i}}(A)$ and $b_{i}=F_{X_{i}}(B)$. Since $Y$ is defined as a minimum, the lower bound does not provide an issue; we simply require that each of the $X_{i}$ be greater than $A$. The upper bound requires constraining the space such that it is possible for one of the $X_{i}$’s to exceed $B$ as long as the other does not. As mentioned above, rejection sampling may be used, but bounds resulting in low-density regions suffer from long running times. Our method draws samples from the unrestricted unit square and maps the results onto the restricted subspace.

The method is composed of the following steps. First generate a sample $(u_{1},u_{2})$ from the uniform distribution over the unit square. Then, based on the value of $u_{2}$, determine which region $R_{s}\subset C$ this point falls into and transform the coordinates using $u_{j}^{\prime}=g_{s}(u_{j})$. These new coordinates fall within $C^{\prime}$ with the same probability distribution that rejection sampling would give but satisfies the constraints with only a single draw needed. After the draw, apply the inverse CDFs to the $(u_{1}^{\prime},u_{2}^{\prime})$ to transform to the original two random variables with $X_{j}=F_{j}^{-1}(u_{j}^{\prime})$. Lastly, take the minimum to obtain a sample from our desired distribution: $Y=\min(X_{1},X_{2})$. Figure 1 shows a diagram of the two spaces and their divisions for some example values, graying out the regions where samples in $C^{\prime}$ cannot be located.

To partition the restricted space $C^{\prime}$ into regions, first note that there are three ways which the $X_{i}$ could satisfy the bounds: only $X_{1}$ is less than $B$, both $X_{1}$ and $X_{2}$ are less than $B$, or only $X_{2}$ is less than $B$. $B$ is the only bound mentioned because in all scenarios both random variables must be greater than $A$. Create an ordered list of the index combinations for $X_{i}$ which are less than $B$ in each possibility: $S=\{\{1\},\{1,2\},\{2\}\}$. Assign coordinates $(u_{1}^{\prime},u_{2}^{\prime})$ to the axes of our unit square, and divide it into three regions $R_{s}^{\prime}$ for each $s\in S$ where

[TABLE]

This defines the region $R_{s}^{\prime}$ to be the one where the coordinates with indices in $s\in S$ lie within their boundaries, and the other coordinates exceed them. We can see that this division of the square forms a nonoverlapping cover of $C^{\prime}$, i.e. $\bigcap\limits_{s\in S}R_{s}^{\prime}=\emptyset$ and $\bigcup\limits_{s\in S}R_{s}^{\prime}=C^{\prime}$. In order to map into these sets without disrupting the probability distribution, we need to know the relative volumes of these regions. We see that for region $R_{s}^{\prime}$, its volume is $V_{R_{s}^{\prime}}=\prod\limits_{j\in s}(b_{j}-a_{j})\prod\limits_{k\notin s}(1-b_{k})$. Furthermore, the volume of the region we cut out of the unit square to create the region $C^{\prime}$ is $V_{\text{out}}=\prod\limits_{j=1}^{2}(1-b_{j})$, so we know that the total volume of region $C^{\prime}$ is $V_{C^{\prime}}=1-V_{\text{out}}$. Thus, the fraction of the total volume which is occupied by $R_{s}^{\prime}$ is $f_{s}=\frac{V_{R_{s}^{\prime}}}{V_{C^{\prime}}}$. Because we are working with a uniform distribution, we need this fraction to be the probability that our final uniform sample came from region $R_{s}^{\prime}$.

We want to start with a uniform draw over an unrestricted unit square which we will call $C$. Call the coordinates of this sample $(u_{1},u_{2})$, where $\forall i\in\{1,2\},0<u_{i}<1$. We need to find a mapping from $C$ to $C^{\prime}$ which preserves the relative probabilities of falling into the regions $R_{s}^{\prime}$. To do this, we can cut along coordinate $u_{2}$ to create three regions $R_{s}$ with volumes $f_{s}$ in our new space so that we can map region $R_{s}$ into $R_{s}^{\prime}$ and preserve this probability. Thus, in our unrestricted space, we can define:

[TABLE]

so that the total volume of $R_{s}$ is $f_{s}$. Once we have this division, it is easy to create the mapping from $R_{s}$ to $R_{s}^{\prime}$: $g_{s}(u_{j})=u_{j}^{\prime}$ where

[TABLE]

where the function $g_{s}(u_{j})$ indicates that this $g_{s}$ is only to be used to map between points in $R_{s}$ and $R_{s}^{\prime}$, where $(u_{1},u_{2})\in R_{s}$ if $\sum\limits_{j=1}^{s-1}f_{j}<u_{2}<\sum\limits_{j=1}^{s}f_{j}$. Define the final mapping to be $g(u_{1},u_{2})=(g_{s}(u_{1}),g_{s}(u_{2}))$ when $(u_{1},u_{2})\in R_{s}$.

3 General Procedure

The following is the general procedure for $N$ random variables of which we want the truncated $k$-th order statistic. Denote the unrestricted $N$-dimensional unit hypercube from which we will draw our original uniform sample as $C=[0,1]^{N}$, and let points in this space have coordinates given by $(u_{1},\cdots,u_{N})$. Let the bound restricted hypercube be $C^{\prime}$ with points having coordinates $(u_{1}^{\prime},\cdots,u_{N}^{\prime})$. Our original bounds, $A$ and $B$, apply to the final $Y$ and not the uniform samples. By applying the CDFs, $a_{i}=F_{X_{i}}(A)$ and $b_{i}=F_{X_{i}}(B)$, we obtain the appropriate limits on $u_{i}^{\prime}$. Using these definitions, we can see that the restricted region $C^{\prime}$ from which we need to generate uniform samples can be written as:

[TABLE]

Our procedure divides $C$ and $C^{\prime}$ into the corresponding regions, $R_{s}$ and $R_{s}^{\prime}$, and determines the mapping of points from $R_{s}$ to $R_{s}^{\prime}$. After the mapping we can use ordinary inverse transform sampling to generate our draws from $Y$.

3.1 Space Segmentation

Each draw must satisfy the constraint that at least $k$ of the $X_{i}$ must be less $B$ in order for $Y$ to be less than $B$, and at least $N-k+1$ of them must be greater than $A$ in order for $Y$ to be greater than $A$. To accomplish the task, we enumerate all the possible combinations of allowed $X_{i}$ orderings. To accomplish this, let us ease the notation of each $X_{i}$ to be represented by its index $i$, and let $I=\{1,\cdots,N\}$ as the set of all possible indices. Then, we can define two sets composed of subsets of $I$

[TABLE]

where $\mathcal{P}(I)$ is the power set of $I$. The first set, $S_{A}$, is the list of possible combinations of the $X_{i}$ which could fulfill the $A$ bound and likewise for $S_{B}$. However, not every pair of $s_{A}\in S_{A}$ and $s_{B}\in S_{B}$ is a valid combination. For instance, if $N=4,k=3$ and we have two sets $s_{A}=\{1,2\},s_{B}=\{1,2,4\}$, then we see that $3\notin s_{A}$ and $3\notin s_{B}$. This means that $X_{3}<A$ and $X_{3}>B$ but since $A<B$, then we have a contradiction. In order to avoid such problems, we must impose two restrictions on the pair: $s_{A}\cup s_{B}=I$ so that no index is omitted and $s_{A}\cap s_{B}\neq\emptyset$ so that there is a value satisfying the bounds to be our $k$-th order statistic. We define the set of allowed index combinations to be:

[TABLE]

Assume that this set is given an order and can therefore be indexed by integers between $1$ and $|S|$.

Using the above index list, we can divide $C^{\prime}$ into smaller regions $R_{s}^{\prime}$ indexed by the elements $s\in S$, where

[TABLE]

Thus, $R_{s}^{\prime}$ is the region of the unit hypercube in which all $X_{i}$ for $i\in s_{A}$ are greater than $A$ and all $X_{i}$ for $i\in s_{B}$ are less than $B$. It is easy to see that $R_{s_{1}}^{\prime}$ and $R_{s_{2}}^{\prime}$ are non-overlapping for $s_{1}\neq s_{2}$ since changing any index’s presence in either set changes the allowed range of values. Additionally, the union of these regions is the total restricted region $C^{\prime}$ since there are no further combinations which give us a $k$-th order statistic in the allowed range of values.

In defining the regions $R_{s}\subset C$ corresponding to the $R_{s}^{\prime}$, we ensure that the probability distribution after the mapping is uniform over $C^{\prime}$ to conduct inverse transform sampling. The probability densities over both $C$ and $C^{\prime}$ follow a multivariate uniform distribution. Because the density function is inversely proportional to the volume, using our mapping we can maintain the relative volumes of the regions being mapped to one another. Then the volume of the regions in $C$, $V_{R_{s}}$, are equal to the normalized volumes in the restricted regions in $C^{\prime}$, $f_{s}$. The normalizing of the volumes begins with the un-normalized volume of $R_{s}^{\prime}$, $V_{R_{s}^{\prime}}$, which is given by

[TABLE]

Then, the normalized volume is the fraction of the allowed volume occupied by the region $R_{s}^{\prime}$ and is given by $f_{s}=\frac{V_{R_{s}^{\prime}}}{\sum\limits_{\sigma\in S}V_{R_{\sigma}^{\prime}}}$.

Then, we can define the regions $R_{s}$ in $C$ so that they have volume $f_{s}$

[TABLE]

Clearly, our region $R_{s}$ has volume $f_{s}$, and, similarly to $R_{s}^{\prime}$, $R_{s_{1}}$ and $R_{s_{2}}$ are non-overlapping and the union of all of the $R_{s}$ comprises the whole space $C$.

3.2 Mapping Between Spaces

Now we have defined a set of allowed orderings of the random variables as $S$ whose elements $s$ identify regions $R_{s}\subset C$ and $R_{s}^{\prime}\subset C^{\prime}$. Once we have our uniformly sampled point $(u_{1},\cdots,u_{N})$ from $C$, we use the above definition of the regions to determine which $R_{s}$ our point falls into. Once we know this $s$, we can use the following mapping to translate the point into $C^{\prime}$ coordinates $(u_{1}^{\prime},\cdots,u_{N}^{\prime})$. For coordinate index $j$ and region index $s=(s_{A},s_{B})$, we have the mapping

[TABLE]

After performing this mapping on each of the $N$ coordinates in our original sampled point, we obtain $(u_{1}^{\prime},\cdots,u_{N}^{\prime})$. Using the inverse CDFs, we can use these coordinates to obtain samples from the $X_{i}$: $(F_{X_{1}}^{-1}(u_{1}^{\prime}),\cdots,F_{X_{N}}^{-1}(u_{N}^{\prime}))$. Then our final answer is the $k$-th largest of these values, which due to our procedure is guaranteed to be between $A$ and $B$ without biasing the distribution.

4 Simulation

We can demonstrate the validity and computational speed of the above method by simulating draws in R.111R code for performing these procedures is available upon request. For the following simulations, we choose $N=5$ and $k=3$, and we use the following probability distributions for the individual $X_{i}$:

• •

$X_{1}\sim\operatorname{Cauchy}(x_{0}=5,\gamma=1)$

• •

$X_{2}\sim\operatorname{Normal}(\mu=6,\sigma=2)$

• •

$X_{3}\sim\operatorname{Logistic}(\mu=3,s=2)$

• •

$X_{4}\sim\operatorname{Weibull}(\lambda=10,k=1.5)$

• •

$X_{5}\sim\operatorname{Uniform}(a=-5,b=20)$

We perform draws using both our new method and ordinary rejection sampling. First, we can compare the empirical CDFs generated from $10000$ draws using both methods on a variety of ranges.

In Figure 2, we see that there is great agreement between the CDFs as determined by the two methods, to the level expected for $10000$ draws. Once certain that our method is accurately drawing from the truncated order statistic distribution, we can assess how long it takes relative to rejection sampling. In determining timing, we used a Dell Latitude E7450 running Windows 10 Pro with a 2.6 GHz Intel Core i7 processor and 16GB 1600 MHz DDR3 memory. For the timing comparison, we used six different intervals as bounds and sampled from each interval $1$, $10$, $100$, $1000$, and $10000$ times, comparing the time it takes to draw all of the samples.

In Figure 3(a), we show that both rejection sampling and our new method behave similarly as the number of samples drawn is increased. This behavior is approximately linear, which makes sense due to the fact that any overhead time would be negligible in comparison with the amount of time per draw, though this overhead time for our new method can be seen in the relatively higher time taken for only one sample. In addition to the number of samples drawn, we need to consider the probability density between the bounds being considered.

Figure 3(b) shows how the execution time varies as the bounds are changed. The bounds themselves are not of interest here since an interval of fixed length can be very easy to sample from if near the peak of the distribution or hard to do if in the tail. Instead, the area of the PDF in that interval is what controls how long rejection sampling may take. The new method is less efficient than rejection sampling when the PDF area is above 5-10%, but below that it can take considerably less time. The extra time taken by the new method is largely due to the extra computations involved in the mapping and determining the mapping groups. The running time of the new method should not vary with the PDF area since every draw leads to a valid sample, which is supported by the figure, but rejection sampling varies by multiple orders of magnitude as the area gets smaller due to the number of draws needed to find a valid one. If a fraction $f$ of the total PDF area is available within the bounds, then it will take $\frac{1}{f}$ draws to find a valid sample on average from rejection sampling. The smallest PDF areas used here were  0.1% of the total area, but as this gets lower (which could occur for narrow ranges or near the tails of a distribution) the time taken should continue to increase dramatically.

5 Conclusion

By mapping unrestricted uniform draws from the $N$-dimensional unit hypercube onto a restricted subspace based on given lower and upper bounds, we can utilize inverse transform sampling to simulate draws from the $k$-th order statistic of $N$ random variables with known CDFs and inverse CDFs. This procedure results in an unbiased draw from this distribution, and in the case where the PDF density in the region of interest is low, the speed of execution is much faster than rejection sampling. This method could be used in place of rejection sampling to ensure that the execution will take a known amount of time if the density of the PDF is unknown in the bounded volume. A practical example would be in drawing from the tail of the minimum or maximum of a set of random variables which are not identically distributed; in such a case, this method would dramatically outperform rejection sampling.