Simple rejection Monte Carlo algorithm and its application to multivariate statistical inference

Statistical inference plays a crucial role in all scientific fields, serving to provide theoretical underpinnings for quantitative models of uncertainty and validate data-based models. Modern statistical inference methods based on sampling enable us to comprehend complex phenomena (Warne et al., 2018). Sampling can be categorized into two types: population sampling and data set sampling. The former is often referred to as Monte Carlo or random simulation, which can be further divided into distribution function-based sampling (Meng et al., 2022) and density function-based sampling (Chen and Yang, 2006; Andrieu and Freitas, 2003). The latter is commonly known as subsampling (Ai et al., 2021; Li et al., 2022), which also falls under the category of Monte Carlo when we consider a data set as a special representation of the population.

Monte Carlo sampling is commonly employed for generating random numbers of random variables or random vectors, which are then utilized to estimate unknown mathematical expectations and their confidence intervals. This well-known algorithm, known as the Monte Carlo algorithm, can also provide an estimation of the complete distribution (L’Ecuyer et al., 2022). The fundamental concept behind the Monte Carlo algorithm lies in utilizing sampled data to perform various calculations and feature analyses. Generating random numbers that adhere to a probability distribution model becomes a crucial technology for these algorithms, often accomplished by computers. There exist numerous options for sampling methods; among them, Markov Chain Monte Carlo (MCMC) (Robert and Casella, 2010) is widely adopted due to its ease of implementation. MCMC technique serves as the gold standard technique for Bayesian inference (Nemeth and Fearnhead, 2021), frequently used to extract samples from challenging posterior distributions (Vyner et al., 2023). However, correlations exist between data generated by MCMC, posing challenges in convergence assessment. In recent years, the Rejection Monte Carlo (RMC) algorithm has emerged sporadically in literature under alternative names such as screening sampling or reject-accept sampling algorithms (Warne et al., 2018; Martino and M´ıguez, 2011; Ma et al., 2022).

Compared to MCMC technology, there is a relatively limited amount of literature on RMC technology. RMC is commonly employed for constructing complex Monte Carlo algorithms, but it itself serves as a standard Monte Carlo technique that filters samples generated by the proposal distribution through the density function to obtain samples from the target distribution (Chen and Yang, 2006; Maryino and M´ıguez, 2011). In comparison with MCMC, RMC possesses unique advantages (Warne et al., 2018): (1) RMC is user-friendly; (2) The generated samples are independent and equally distributed; (3) Calculation efficiency remains unaffected by calculation parameters, particularly without requiring heuristic definition of transition probability. However, the drawbacks of RMC are also evident as it can be challenging to determine an appropriate recommendation sampling distribution. When the recommendation distribution significantly deviates from the target distribution, computational efficiency becomes low. Therefore, previous studies have primarily focused on sampling efficiency and determining efficient recommendation distributions (Chen and Yang, 2006; Maryino and M´ıguez, 2011). Nevertheless, with continuous advancements in computer technology, efficiency no longer poses a key constraint in many cases. The flexibility of RMC technology towards suggestion distributions has been further enhanced thereby making its application more widespread. Consequently,this paper provides enlightening insights into the convenient learning perspectives and program implementation aspects of RMC techniques to enable arbitrary model sampling using RMC methods—especially for multivariate models—and explores their versatility and practical value.

The RMC algorithms are a versatile approach for generating independent samples from a target density (Martino and M´ıguez, 2011). Let $f(x)$ be the probability density function of a target random variable $X$, and let $f_{0}(x)$ the suggested density function that can be easily sampled, with the same support supp(f) as the function $f(x)$. Additionally, let $c$ a constant satisfying the following condition:

The RMC algorithms are typically delineated in the literature as follows.

A merit function, such as $cf_{0}(x)$ in Inequality $(1)$ plays a crucial role in the efficiency of the sampling. When the deviation between the proposed model ${f_{0}}(x)$ and $f(x)$ is significant and $c$ is chosen excessively large, the difference between the resulting merit function and objective function, $cf_{0}(x)-f(x)$, becomes substantial. This leads to a low acceptance probability $P\{y\leqslant\frac{{f(x)}}{{c{f_{0}}(x)}}\}$, causing most points to be rejected during random selection and ultimately reducing sampling efficiency.

The construction of $f_{0}(x)$ is typically designed to enhance sampling efficiency by closely approximating $f(x)$, while the constant $c$ is determined as the $c=\arg{\min_{\text{d}}}\{d{f_{0}}(x)\geqslant f(x),x\in{\text{supp}}(f)\}$.

Theorem 1 Let the value $x$ of random variable X be a random number generated by the Algorithm $1$, then X $\sim f(x)$.

Proof The density functions of random numbers $x$ and $y$ are denoted as ${f_{0}}(x)$ and ${f_{1}}(y)=\left\{{\begin{array}[]{*{20}{c}}{1,}&{y\in[0,1]}\\ {0,}&{\rm other.}\end{array}}\right.$, respectively, according to the Algorithm 1. As the processes generating $x$ and $y$ are mutually independent, the joint probability density function of $x$ and $y$ can be expressed as follows:

The cumulative probability function of random number $x$:

The identification of the optimal merit function constitutes a fundamental aspect in existing RMC literature (Chen and Yang, 2006; Martino and Mígues, 2011). However, its practical implementation poses significant challenges, thereby becoming a bottleneck that hinders the progress of the RMC algorithm. According to Algorithm 1, the construction of the proposed distribution is unrestricted and may deviate from the target distribution. This discrepancy inevitably leads to inefficient sampling; nevertheless, advancements in computer technology have compensated for this limitation. Consequently, contemporary approaches favor either uniform distribution models or segmented uniform distribution models as suitable candidates for the proposed distributions. By employing uniformly distributed or segmented uniformly distributed merit functions, it is possible to significantly enhance both versatility and comprehensibility within RMC techniques while effortlessly facilitating computer-based sampling across arbitrary distribution models.

An intuitive understanding of sampling from the target density is to allocate more samples in regions with higher density values and fewer samples in regions with lower density values, proportionally matching the value of the density function. This can be easily achieved by employing a uniform distribution. For instance, when considering sampling from a given density function

The uniform distribution can be employed to evenly distribute points across the rectangular area $[\pi/4,3\pi/4]\times[0,1.1]$. Subsequently, the points located between the density function and the x-axis effectively reflect the characteristics of the density function and are thus considered valid. As depicted in Figure 1, these accepted points correspond to samples based on their horizontal coordinates.

The previous analysis suggests using a simple technique called Simple Rejection Monte Carlo (SRMC) to sample from a model $f(x)$, which can be described as Algorithm 2.

Given the ease of drawing samples uniformly from the support set ${\rm suff}(f)$ of function $f(x)$, Algorithm 2 can be straightforwardly extended to a Simple Rejection Monte Carlo (SRMC) technique for sampling from a multivariate random variable model $f({x_{1}},{x_{2}},\cdots,{x_{d}})$, as described in algorithm 3.

As computing technology continues to advance, computational efficiency is often no longer the primary constraint, and the implementability of algorithms assumes greater significance. Algorithm 2 and Algorithm 3 can be easily implemented programmatically, thereby facilitating intelligent computation. The evaluation of algorithm strengths and weaknesses can be measured in terms of accuracy and program runtime. The following specific examples demonstrate that RMC out-sampling for multivariate distributions can be readily implemented, while considering bias and program runtime as indicators of accuracy and feasibility. All programs are implemented using the R language (Meng, 2022).

Example 1 Let the density function of a two-dimensional Gaussian random vector $(X,Y)$ be

where $\mu=\left({\begin{array}[]{*{20}{c}}x\\ y\end{array}}\right)\in{R^{2}}$, $\sum=\left({\begin{array}[]{*{20}{c}}1&{0.2}\\ {0.2}&1\end{array}}\right)$. Please generate a sample of the random vector by RMC technology.

According to Algorithm 3, the square region $[-5,5]\times[-5,5]$ is selected with $c=\max f(x,y)=0.1657$. The scatter plot of the sampling is presented in Figure 2. The sample sizes from left to right are: 1000, 10000 and 100000; corresponding program running times are 0.59, 5.19 and 51.65 seconds respectively. The correlation coefficients calculated from the extracted samples are found to be 0.237, 0.189 and 0.202 respectively, which exhibit slight deviation from the theoretical correlation coefficient of 0.2; moreover, as the number of samples increases, the absolute value of this deviation decreases accordingly.The scatterplot (Figure 2) visually demonstrates the structural relationship between variables where most values fall within the square region $[-4,4]\times[-4,4]$, while occurrences outside this range can be considered negligible.

Probability calculations play a crucial role in statistical inference, often necessitating the use of integral computations. However, solving these integrals analytically can be challenging, thus making sampling methods an invaluable tool for their resolution. This is exemplified below through the application of ordinary dual integration.

Example 2 Calculate the dual integral $\iint_{{\rm D}}{xy{{\rm d}}x{{\rm d}}y}$, where the region D is bounded by ${y^{2}}=x$, $x$ axis and $y=x-2$, and $D=\{(x,y)\mid y\geq 0\}.$

Construct a rectangular area $S=[0,4]\times[0,2]\supseteq D$. Then the integrated function $xy$ is easily integrated over $S$, and the function

can be viewed as a density function of a random vector $(X,Y)$. Then the required dual integral can be transformed into

where ${\rm I}(x,y)$ is a representative function of the set D, i.e., when $(x,y)\in D$, ${\rm I}(x,y)=1$ ; otherwise, ${\rm I}(x,y)=0$. Since $f(x,y)$ is a density function on S, it follows that$\iint_{S}I(x,y)f(x,y)dxdy=E(I(X,Y))\approx\frac{{\#}\{(x,y)\mid(x,y)\in D\}}{{\#\{(x,y)\mid(x,y)\in S\}}},$ where the symbol $\#$ denotes the number of elements in the corresponding set. If one considers a uniform distribution over the region S, then $\iint_{S}{xy{{\rm d}}x{{\rm d}}y}=8\iint_{S}xy\frac{1}{8}{\rm d}x{\rm d}y=8E(XY)\approx 8\frac{1}{N}\sum_{i=1}^{N}x_{i}y_{i}.$ Therefor, the target integral can be approximated as follows.

From Equation (5) and Algorithm 3, in order to complete the computation of the objective integral, first a uniform sampling should be performed on the region S, and then the resulting uniform samples should be screened by the density function $f(x,y)$. Assume that the set of uniform random numbers drawn first is denoted as S1 and the set of random numbers obtained by subsequent screening is S2. Then the objective integral can be completed by calculating $\frac{1}{N}\sum\limits_{i=1}^{N}{{x_{i}}{y_{i}}}$ and $\frac{{\#}\{(x,y)\mid(x,y)\in D\}}{{\#\{(x,y)\mid(x,y)\in S\}}}$ from S1 and S2, respectively. The results were calculated using R software programming as follows, with 10 runs for each sample size, and the results are averaged.

Previous literature primarily focuses on rejection sampling algorithms for one-dimensional random variables, which involve constructing a proposed sampling model and optimal function that are both easy to sample from and closely approximate the target model. However, despite its high sampling efficiency, the constructed proposed model lacks operability, limiting its intelligent implementation in statistical inference. This paper highlights the versatility and ease of operation provided by the RMC technique and introduces the SRMC technique as a solution. The proposed technique is applicable to both one-dimensional and multivariate random variables, offering generalizability and comprehensibility. Specific examples demonstrate that this algorithm can be easily programmed and implemented in multivariate statistical analysis with high accuracy. Consequently, it will undoubtedly facilitate widespread adoption of the SRMC technique while providing an alternative method for intelligent implementation of multivariate statistical inference.