{'fraction_non_alphanumeric': 0.03688866633283258, 'fraction_numerical': 0.009236076336727314, 'mean_word_length': 4.272367263127017, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Fast Simulation Method for the Sum of Subexponential Distributions

Estimating the probability that a sum of random variables (RVs) exceeds a given threshold is a well-known challenging problem. Closed-form expression of the sum distribution is usually intractable and presents an open problem. A crude Monte Carlo (MC) simulation is the standard technique for the estimation of this type of probability. However, this approach is computationally expensive especially when dealing with rare events (i.e events with very small probabilities). Importance Sampling (IS) is an alternative approach which effectively improves the computational efficiency of the MC simulation. In this paper, we develop a general framework based on IS approach for the efficient estimation of the probability that the sum of independent and not necessarily identically distributed heavy-tailed RVs exceeds a given threshold. The proposed IS approach is based on constructing a new sampling distribution by twisting the hazard rate of the original underlying distribution of each component in the summation. A minmax approach is carried out for the determination of the twisting parameter, for any given threshold. Moreover, using this minmax optimal choice, the estimation of the probability of interest is shown to be asymptotically optimal as the threshold goes to infinity. We also offer some selected simulation results illustrating first the efficiency of the proposed IS approach compared to the naive MC simulation. The near-optimality of the minmax approach is then numerically analyzed.


I. INTRODUCTION
The performance analysis of communication systems is generally associated with the investigation of the statistics of sums of Random Variables (RVs). For instance, when diversity techniques such as Maximum ratio Combining (MRC) and Equal Gain combining (EGC) are performed, the resulting received signal-to-noise-ratio (SNR) is modeled by a sum of fading variates [1].
Unfortunately, the statistics of the sum distribution for most of the challenging problems are generally intractable and unknown. Monte Carlo (MC) simulation is the standard technique to estimate the probability that a sum of RVs exceeds a given threshold. However, this approach requires an extensive computational work to estimate extremely small probabilities. Importance Sampling (IS) is an alternative approach which aims to improve the computational efficiency of the naive MC simulation technique [2]. The basic idea behind IS technique is to change the underlying sampling distribution in a way to achieve a substantial variance reduction of the IS estimator. Many research efforts have been carried out to propose efficient IS algorithms. For instance, among the first works in the digital communication field, the authors in [3] and [4] proposed methods based respectively on scaling the variance and shifting the mean of the original probability measure. An extension of [3] was performed in [5] where a composite IS technique was derived. In [6], the asymptotic efficiency of five different IS techniques was studied for the estimation of the Bit Error Rate (BER) in digital communication systems with Gaussian input. Exponential twisting, derived from the large deviation theory, is an interesting IS change of measure technique since in most of the cases it yields "optimal" asymptotic results [7] [8]. For instance, this technique was used to estimate the BER of direct-detection optical systems employing avalanche photodiode (APD) receivers in [9].
The exponential twisting change of measure is feasible only with distributions having finite Moment Generating Function (MGF). Thus, in the heavy-tailed setting where the MGF is infinite, it is not possible to use the exponential twisting method. However, many heavy-tailed distributions, such as the Log-normal and the Weibull (with shape parameter less than 1) RVs, are frequently encountered in various applications. In cellular mobile communication systems, the Co-Channel Interference (CCI) power which arises due for instanceto the neighboring cells that use the same frequency is generally modeled as a sum of Log-normal (SLN) RVs [10]. Besides, the Log-normal distribution is also used to model the largescale fading in the ultra-wideband (UWB) communications [11], and the weak-to-moderate turbulence channels in freespace optical communication channels [12]. Recently, the Weibull fading has also received an increasing attention since it exhibits a good fit to experimental fading data for both indoor and outdoor environment [13], [14], [15]. Various closedform approximations of the sum of Log-normal RVs [16] [17] [18] [19] and the sum of Weibull RVs [20] [21] [22] have been extensively developed. These approximations are not generic and depend on the problem under consideration. Hence, a lot of research efforts have focused in developing a generic efficient IS technique dealing with distributions in the heavy-tailed class. In [23], two efficient techniques for the estimation of the probability that the sum of subexponential RVs exceeds a given threshold have been presented. The first one is based on conditional MC, whereas the second is based on considering a new probability measure which is heavier than the underlying distribution. In [24], a transform likelihood ratio approach was derived to switch the heavy-tailed problem into an equivalent light-tailed one. The authors in [25] have developed an efficient fast simulation method for estimating a sum of independent and identically distributed (i.i.d) RVs with subexponential decay. Their approach is based on twisting the hazard rate of the original probability measure of each component in the summation.
In this paper, inspired by [25], we develop a general approach based on hazard rate twisting to efficiently estimate the probability that a sum of independent and non-identically distributed heavy-tailed RVs exceeds a given threshold. The twisting parameter is determined through a minmax approach which first ensures a nearly optimal computational gain in terms of the number of simulation runs and second leads to an asymptotic optimality criterion. The rest of the paper is organized as follows. In section II, we state the problem setting and enumerate the main contributions. In Section III, a minmax hazard rate twisting approach is introduced with an emphasis on the general procedure leading to an efficient choice of the twisting parameter. Moreover, the asymptotic optimality criterion using this proposed IS approach is verified. In Section IV, two applications of distributions with subexponential decay are studied. In Section V, a substantial computational gain of the proposed IS technique is analyzed and shown through various selected simulation results.

II. MATHEMATICAL BACKGROUND A. Problem Setting
Let X 1 , X 2 , ..., X N be a sequence of independent but not necessarily identically distributed positive RVs. Let us denote the Probability Density Function (PDF) of each X i by f i (x), i = 1, 2, ..., N . Our objective is to efficiently estimate for a sufficiently large threshold γ th . We focus on heavy-tailed distributions, i.e distributions which exhibit slower decays than the exponential distribution. Formally, a distribution of a RV X is said to be heavy-tailed if In practice, all commonly used heavy-tailed distributions belong to the subclass of subexponential distributions. In fact, a distribution of a RV X is said to be subexponential if where F (x) is the Complementary Cumulative Distribution Function (CCDF) of X, and F * n (x) is the CCDF of the sum of n i.i.d RVs with distribution F . Examples of such subexponential distributions are: the Log-normal distribution, and the Weibull distribution with shape parameter less than 1. The readers are referred to [25] for more discussion about subexponential distributions.
The standard technique to estimate α is to use the naive MC estimator defined aŝ where M is the number of simulation runs, and 1 (·) defines the indicator function. It is widely known that the naive MC simulation is extensively expensive for the estimation of rare events. In fact, from the Central Limit Theorem (CLT), it can be shown that the MC estimation with 10% relative precision requires more than 100/α simulation runs. Hence the number of samples to estimate a probability of order 10 −9 should be more than 10 11 , with an accuracy requirement of 10%.
Consequently, there is a detrimental need to improve the computational efficiency of the MC simulation.
B. Importance Sampling IS is a variance reduction technique which aims to increase the computational efficiency of the naive MC simulation [2]. The general concept of IS is to construct an unbiased estimator of the desired probability with much smaller variance than the naive estimator. In fact, this technique is based on performing a suitable change of the sampling distribution as follows where the expectation is taken with respect to the new probability measure p * under which the PDF of each X i is g i , and L is the likelihood ratio defined as The idea behind this change of measure is to enhance sampling important points which have more impact on the desired probability. Hence, emphasizing that important points are sampled frequently will result in a decrease of the variance of the IS estimator. The new IS estimator is defined aŝ Generally, it is not obvious how to construct a new probability measure which results in decreasing the variance of the IS estimator and hence increasing the computational efficiency. Besides, it is necessary to define some performance metrics which measure the goodness and the pertinence of the IS estimator. Bounded relative error, asymptotic optimality, and bounded likelihood ratio are useful indicators to characterize a good change of probability measure [2]. Generally, it is difficult to achieve the bounded relative error criterion, whereas the asymptotic optimality could be shown if one choose an appropriate probability measure g i . Let us consider the sequence of the RVs {T γ th } defined as From the non-negativity of the variance of T γ th , we get Applying the logarithm on both side, we conclude that for all p * we have Hence, we say that α is asymptotically optimally estimated under the probability measure p * if the above equation holds with equality as γ th → +∞, that is It is important to note that the naive simulation is not asymptotically optimal for the estimation of α since the ratio in (11) is equal to 1.
The exponential twisting technique, which is derived from the large deviation theory, is the main IS framework dealing with light-tailed distributions, that is distributions whose tails decay at an exponential rate or faster. The exponential twisting by an amount θ > 0 is given by where M Xi (θ) denotes the moment generating function (MGF) of the RV X i . In most of the cases, this technique achieves the asymptotic optimality criterion given in (11) [8].
In the heavy-tailed setting, the exponential twisting change of measure is not feasible and alternative techniques are needed. In fact, the MGFs are infinite for distributions with heavy tails. In [25], an efficient IS technique was developed for the estimation of α in the case of i.i.d sum of RVs with subexponential decay. Their idea was based on twisting the hazard rate of each component in the summation S N by a quantity 0 < θ < 1. Let us define the hazard rate λ i (·) associated to the RV X i as where F i (·) is the CDF of X i , i = 1, ..., N . Besides, we define also the hazard function as From (13) and (14), the PDF of X i is related to the hazard rate and function as The change of probability measure is obtained by twisting the hazard rate of the underlying distribution by a quantity 0 < θ < 1 as follows Consequently, the RV T γ th has the following expression For heavy-tailed distributions, the hazard rate twisting based approach plays the same role as the exponential twisting technique in the light-tailed setting. In [26], the authors emphasize the central role played by hazard rate technique for the estimation of small probabilities that a general function containing both light and heavy-tailed distributions exceeds a given threshold. In fact, by developing log-asymptotic expressions for both the probability of interest and the second moment of T γ th , they have proved that α is asymptotically optimally estimated. The equivalence between the hazard rate and the exponential twisting techniques is also emphasized in [27] where a suitable hazard function transformation is used, in the case of a sum of i.i.d subexponential distibutions, to switch from a heavy-tailed problem to a light-tailed one where the exponential twisting could be used.

C. Main Contributions
A primordial question when using either exponential or hazard rate twisting techniques is the choice of the twisting parameter θ. The selection of this parameter should be performed in a way to ensure a maximum reduction of the second moment of T γ th , and hence result in a maximum amount of computational gain. Unfortunately, this is not feasible in general since E θ T 2 γ th ( E θ [·] denotes the expectation under the IS probability measure ) is typically not known in a closed form. The commonly used procedure to determine θ starts by deriving a close upper bound on E θ T 2 γ th and then finding the value of θ which minimizes that upper bound. For the exponential twisting, this upper bound is easily obtained using (12) and (6) Then, the value of θ = θ * selected to minimize the upper bound is satisfying In the hazard rate twisting setting, the determination of θ * is not as straightforward as for the exponential twisting approach. In fact, the upper bound on the second moment is not easy to obtain. In [25], the i.i.d sum of subexponential distributions is considered. The determination of the twisting parameter was done via the derivation of an upper bound on the second moment of T γ th which holds only for a sufficiently large threshold. More precisely, by assuming that the hazard rates are eventually decreasing to zero and are eventually everywhere differentiable, the asymptotic inequality holds for every ǫ > 0 and with N i x i large enough. Then, using the previous asymptotic inequality, an upper bound on E θ T 2 γ th was computed which is minimized when Moreover, they proved in [25] that asymptotic optimality holds by replacing N in (21) by any positive constant. In the present work, we consider a non-trivial generalization of [25] to the case of the sum of independent and non-identically distributed subexponential RVs. Our procedure for the determination of the twisting parameter is performed in two steps. First, we derive an upper bound on the second moment of T γ th through the resolution of a constrained maximization problem on the likelihood ratio. Second, we minimize this upper bound over all possible value of θ which results in the so called minmax optimal twisting parameter θ = θ * . For the class of subexponential distributions, we will see that, under a weaker assumption than the one stated in [25] to derive (20), we are able to characterize the behavior of the solution of the maximization problem and detect the region where the maximum is achieved. In a nutshell, the main contributions of the present paper are: • We develop an optimized hazard rate twisting approach for the estimation of α for the case of the sum of independent and non-identically distributed subexponential RVs. The procedure that we will follow to determine θ is based on a minmax approach. This minmax procedure starts by computing the maximum (the most sharpest upper bound) on the second moment of T γ th for all value of γ th . Then, a simple minimization problem is solved to derive the minmax optimal twisting parameter θ * . Besides, we will see also that this choice of θ is efficient since it almost results in the same computational gain as the unknown optimal value (the value that minimizes the actual second moment of T γ th ). In the particular i.i.d sum, we prove that our minmax twisting parameter is equivalent to the one derived in [25] as γ th goes to infinity. • We prove under some realistic assumptions, which are generally satisfied by distributions with subexponential decays, that α is asymptotically optimally estimated using our minmax approach. • Finally, two applications will be studied to clarify how the procedure is applied, and to validate through numerical results the efficiency of the proposed minmax hazard rate twisting approach. The first application considers the sum of independent Log-normal RVs, and the second one deals with the sum of independent Weibull distributions with shape parameter less than 1. It is important to note that in our approach there is no restriction to consider the sum of a mixture of subexponential distributions belonging to different families.

III. PROPOSED HAZARD RATE TWISTING A. General Approach
Generally, an interesting IS change of probability measure for the estimation of rare events is characterized by the property of uniformly bounded likelihood ratio. This property will result in obtaining an upper bound on the second moment of the RV T γ th . Then, the optimal value of the parameter θ is chosen to minimize that upper bound. More precisely, the procedure of choosing θ is divided into two steps. In the first step, we construct an upper bound on the second moment of E θ (T 2 γ th ) which is achieved by solving the following maximization problem (P): where the likelihood ratio is given as follows Hence, solving the problem (P ) is equivalent to solving the following minimization problem (P ′ ): The resolution of the maximization problem (P ) or equivalently the minimization problem (P ′ ) will be discussed later in the following subsection. By denoting the optimal solution of (P ) by X * 1 , X * 2 , ..., X * N , we have The second step is to minimize (25) to get the optimal twisting parameter θ * . This is a simple minimization problem to solve which results in

B. Asymptotic Optimality Criterion
The value of the twisting given in (26) represents the minmax optimal choice among all values of θ, and for all threshold values. Now, we focus on the asymptotic behavior of the IS estimator as γ th goes to infinity. In particular, we investigate the asymptotic optimality criterion (11) using the twisting parameter θ * given in (26).
The investigation of the asymptotic optimality criterion is based on analyzing the asymptotic behavior of the solution of the minimization problem (P ′ ). Since each hazard function Λ i (·) is an increasing function, it follows that the inequality constraint is satisfied with equality, that is In order to ensure the asymptotic optimality, let us consider the following assumption Assumption 1. For each i ∈ {1, 2, ..., N }, we assume that there exist η i such that the hazard function Λ i (·) is concave in the interval [η i , +∞).
The previous assumption is consistent with all commonly used subexponential distributions such as the Log-normal, the Weibull (with shape parameter less than 1), and the Pareto (with parameter strictly bigger than 1) distributions. In the following lemma, we characterize the behavior of the solution of (P ′ ) for a sufficiently large threshold γ th : Lemma 1. Under Assumption 1, there exists a fixed index i 0 ∈ {1, 2, ..., N } such that the minimizers of (P ′ ) satisfy for a sufficiently large γ th and hence as γ th → +∞, we have Proof: Let us consider S(N, γ th ) the set of all feasible solutions: Through the use of Assumption 1, the objective function of (P ′ ) is concave on the subset: Thus, the minimum of the objective function of (P ′ ) overS(N, γ th ) is achieved in at least one of its extreme points. More precisely, the extreme points ofS(N, γ th ) are e 1 , e 2 , ..., e N such that e i = (η 1 , η 2 , ..., η i−1 , γ th − j =i η j , η i+1 , ..., η N ). Therefore the minimum of (P ′ ) over In both cases, there exists at least one index i ∈ {1, 2..., N } such that X * i ≤ η i . In addition, in order to satisfy the equality constraint N i=1 X * i = γ th for a sufficiently large γ th , there should exist an index j ∈ {1, 2, ..., N } such that X j ≥ η j . In order to prove the result in Lemma 1, we proceed iteratively by dimension reduction. In fact, without loss of generality, we assume that X * N ≤ η N (through an index permutation). It follows that Consequently, we can see that we have reduced the number of optimization variables to be N − 1, while we have kept the same structure of the minimization problem (P ′ ) with γ * th,N −1 = γ th − X * N . Hence the previous procedure could be repeated again. In fact, using the same argument as before, there exists another index i ∈ {1, 2, ..., N − 1} such that X * i ≤ η i . Without loss of generality, we assume that i = N −1 which leads to where γ * th,N −2 = γ th − X * N − X * N −1 . After N − 2 steps, we get with X * i ≤ η i , for i = 3, 4...., N , and γ th,2 = γ th − N i=3 X * i . Thus, we end up with a two dimensional minimization problem. Again, there should exist an index i = 2 ( through a possible permutation ) such that X * 2 ≤ η 2 . Therefore, using the equality constraint N i=1 X * i = γ th , we get Since η i , i = 2, 3, ..., N are independent of γ th , it follows Thus, as γ th goes to infinity, we have It is important to note that in the particular i.i.d case, the index i 0 could be any index in {1, 2, ..., N }, and the minimum is achieved in N different points. A direct consequence of Lemma 1 is presented in the following lemma. Lemma 2. Under Assumption 1, the objective function of (P ′ ) has the following asymptotic behavior Proof: Using Lemma 1 and the fact that Λ i0 (γ th ) tends to infinity as γ th increases, we have The remaining work is to prove that Using the fact that Λ i0 (·) is increasing to infinity and concave for inputs bigger than η i0 , then its derivative which is the hazard rate λ i0 (·) is a decreasing function provided that which leads to (46) and then the proof is concluded. Now, we could state the asymptotic optimality theorem Theorem 1. For a sum of independent RVs with subexponential distributions and under Assumption 1, the quantity of interest α is asymptotically optimally estimated using the hazard rate twisting approach with the minmax optimal twisting parameter θ * given in (26).
Proof: In (25), we have derived an upper bound on the second moment of T γ th as and replacing the optimal twisting parameter θ * given in (26), we have By applying the logarithmic function on both side, it follows On the other hand, using the non-negativity of X i , i ∈ {1, 2, ..., N }, we have Note that for sufficiently large γ th , the left and right-hand sides of (51) are negative. Therefore, Finally, using Lemma 2, we have: (54) Through the use of the non-negativity of the variance, we conclude the proof.

Remark 1.
Under the i.i.d case, Assumption 1 is almost equivalent to the one stated in [25]. They assumed also that the hazard rate is converging to zero whereas in our case this is not needed. The previous observation makes our assumption a bit weaker compared to [25]. In addition, our optimal twisting parameter θ * given in (26) tends to the same value (21) derived in [25], as γ th goes to infinity.

C. Generation of the Twisted Distribution
Generally, hazard rate twisting the original PDF of a RV X does not result in a known distribution. One way to generate realizations of X under f θ (·) could be performed via its CDF F θ (·). In fact, it is known that F −1 θ (U ), where U is uniformly distributed RV over [0, 1], has the same distribution as X under the hazard rate twisted PDF [28]. Let us consider a RV X with an underlying PDF f (·) and CDF F (.). From (16), the PDF f θ (·) associated to X with hazard rate λ(·) and hazard function Λ(·) is Replacing λ(·) and Λ(·) by their definitions, we get By a simple integration, the corresponding CDF is given by Finally, a simple computation leads to an exact expression of the CDF inverse of the RV X under the hazard rate twisting technique where F −1 (·) is the CDF inverse of X under the original PDF f (·). A pseudo-code describing all steps to estimate α by our proposed hazard rate twisting approach is described in Algorithm 1.

Algorithm 1
Optimized hazard rate twisting approach for the estimation of α Inputs: M IS , γ th . Outputs:α IS . Find the optimal value of θ as in (26) by solving the minimization problem (P ′ ). for i = 1, ..., M IS do Generate N independent realizations of the uniform distribution over [0, 1]: (17). end for Compute the IS estimator asα IS = 1

IV. APPLICATIONS
We consider two examples of distributions belonging to the class of subexponential distributions: the Log-normal and the Weibull (with shape parameter less than 1) distributions. We will investigate for these two examples the solution of (P ′ ).

A. Weibull Distribution
In this example, the PDF of X i , i = 1, 2, ..., N is where k i > 0 and β i > 0 denotes respectively the shape and the scale parameters. We focus on the case where the shape parameter is strictly less than 1 since it is known that with this choice the Weibull RV is a subexponential distribution. The hazard rate and function for each X i are as follow Let us now investigate the solution of the minimization problem (P ′ ). We could prove through a simple computation that the objective function of (P ′ ) is concave for k i < 1, i = 1, 2, ..., N and hence Assumption 1 is satisfied. In fact, the Hessian H of the objective function at any point which is strictly negative for k i < 1, i = 1, 2, ..., N .
In particular, the objective function is also concave on the convex set S(N, γ th ) = {X = (X 1 , X 2 , ..., X N ) ∈ (R + ) N , such that N i=1 X i = γ th }. Therefore, the solution of (P ′ ) is obtained in one of the extreme points of S(N, γ th ). In other words, the minimum is achieved when where i 0 satisfying It is important to note that for large values of γ th , the index i 0 depends only on the shape and scale parameters and independent of γ th . More precisely, for γ th large enough, it is characterized by Moreover, if there are more than one RV with minimum shape parameter, the index i 0 corresponds to the one with maximum scale parameter. Remark 2. We have described in the previous section a method based on the inverse CDF F −1 θ (·) to generate samples of a RV X under the twisted PDF f θ (·). For the particular Weibull distribution with parameters k and β, the PDF f θ (·) is simply another Weibull distribution with the same shape parameter k and a different scale parameter β ′ as follows where

B. Log-Normal Distribution
The PDF of each X i , i = 1, 2, .., N is given by where µ i and σ i are the mean and the standard deviation of the associated Gaussian RV Y i = log(X i ). In communication, the decibel unit is generally used. Hence, it is more convenient to define a Gaussian RV as Z i = 10 log 10 (X i ) with mean µ i,dB and standard deviation σ i,dB . The relation between the two Gaussian RVs Y i and Z i , i = 1, 2, ..., N , are where ξ = log(10)/10. The expressions of λ i (·) and Λ i (·) are given by where φ(·) and Φ(·) are respectively the PDF and the CDF of a standard Gaussian distribution. In this example, the solution of (P ′ ) is not straightforwardly computed as the Weibull distribution. The difficulty to find out the analytic solution of the minimization problem (P ′ ) arises from the fact that the hazard function for a Log-normal RV does not have the concavity property as for the Weibull distribution.
However, it is known that the hazard function of the Lognormal distribution has the property stated in Assumption 1. Therefore, the minimizers of (P ′ ) satisfies Lemma 1 which states that there exists an index i 0 such that X * i0 is close to γ th whereas the other components are bounded. Hence, since the hazard function Λ i is an increasing function, the index i 0 satisfies for a sufficiently large γ th Thus, for γ th large enough, the index i 0 is independent of γ th and corresponds to Moreover, if there exists another index with a maximum standard deviation, i 0 corresponds to the RV with a maximum mean.

V. SIMULATION RESULTS
In this section, some selected simulation results are shown to compare the naive MC simulation and the proposed IS simulation technique. Two performance metrics will be used to compare these two approaches. The relative error of the naive MC estimator is defined through the use of the CLT as and the relative error of the IS MC estimator is given by where C is the confidence constant equal to 1.96 (for 95% confidence interval), and M MC and M IS are the number of samples for the naive MC and the IS MC simulations, respectively. Note that the use ofα IS in (73) instead of α MC gives a more accurate estimate of the standard deviation ofα MC . For a fixed relative error, we define the efficiency indicator of the IS MC technique compared to the naive MC simulation as The more the efficiency k is large, the more we need samples in the naive MC simulation to reach the relative accuracy given by IS. In other words, the bigger is k, the more efficient is the proposed IS technique.

A. Frequency of Occurrence
As it was mentioned before, a key characteristic of a good IS technique is to emphasize the sampling of important points, i.e the number of realizations satisfying S N ≥ γ th . We define the frequency of occurrence as the number of samples which satisfy S N ≥ γ th . In our first simulation results, we consider the sum of two i.i.d. Log-normal RVs with mean µ dB = 0 dB and standard deviation σ dB = 6 dB. In Table  I,    the IS change of measure where the frequency of realizations belonging to the rare set S N ≥ γ th is almost constant as we increase the threshold. On the other hand, the failure of sampling under the original SLN distribution is clear through its inability to construct realizations in the rare sets. In Table  II, we show the same computation using the sum of two i.i.d Weibull distribution with shape parameter k = 0.5 and scale parameter β = 1. Again, important points are sampled more frequently using the IS technique and their frequencies remains almost constant as we increase the threshold.
To illustrate this statement, we plotted in Fig. 1 the twisted against the original SLN distributions for a fixed threshold γ th = 20. Clearly, we see that twisting the hazard rate of each component in the sum leads to a more heavier twisted PDF. As a consequence, the events which exceed the given threshold are more likely to occur under the twisted PDF than under the original one.

B. Efficiency of the Proposed IS Algorithm
In Fig. 2, the CCDF of the sum of two i.i.d Log-normal RVs is presented using both the naive MC simulation and our IS simulation technique. The inefficiency of the naive simulation is clear in Fig. 2. In fact, a remarkable oscillatory behavior of the naive MC technique is observed using a number of samples M MC = 10 6 for γ th ≥ 25 dB. Besides, as we increase the threshold, the naive MC estimator is almost zero. Indeed, more samples are required in order to overcome this failure and to get a good approximation of the CCDF. The naive technique with M MC = 10 8 is also presented in Fig. 2 and is compared to IS simulation. We point out that both methods coincide and we have a good approximation of the CCDF up to a probability of order 10 −6 . Then, an oscillation of the tail of the CCDF using the naive MC approach is observed, whereas IS technique gives a smooth curve. Thus, our IS technique gives a more accurate result using a less number of samples 5×10 4 , in contrast with 10 8 samples used in the naive simulation. In order to confirm the previous statement, we need to analyze the relative error given by both techniques.  In Fig. 3, we plotted the relative error of the naive and the IS simulations as function of the threshold. We point out a slow variation of the relative error of the naive MC simulation for γ th < 25, then a very rapid increase is observed as we increase the threshold. In fact, in the first region the number of samples is sufficient to guarantee an accurate approximation, whereas in the second region the naive simulation fails to well estimate the CCDF and hence substantial samples are required to ensure a good accuracy, i.e much more than 10 8 realizations. On the other hand, IS technique shows an interesting result in Fig.  3 where the variation of its relative error is extremely slow compared to the naive simulation. Consequently, with M IS much smaller than M MC , our IS approach approximates the CCDF more efficiently than the naive simulation.
In Fig. 4, we plotted the efficiency indicator k as function of the threshold. From this figure, we deduce that the efficiency is increasing rapidly, almost exponentially. Hence, the more we increase the threshold the more efficient is our IS technique.  This result is expected since k is proportional to the number of samples M MC that we need to generate in order to absorb the rapid increase of the relative error of the naive MC simulation, i.e reach the relative accuracy given by the IS approach. Besides, Fig. 4 illustrates also that the IS technique is more efficient for the considered range of probability, i.e k always bigger than 1.
In the second simulation results, we consider the sum of two independent Weibull distribution with same scale parameter β = 1, and with different shape parameters k 1 = 0.4, and k 2 = 0.8. In Fig. 5, Fig. 6, and Fig. 7, we plotted the CCDF, the relative error, and the efficiency, respectively. We note that in this case also, the proposed IS technique gives an accurate and efficient approximation of the CCDF and results in a substantial computational gain.

C. Near-Optimality of the Minmax Twisting Parameter
In our next simulation results, we aim to analyze the sensibility of the second moment of the RV T γ th with respect to the twisting parameter θ. Since our twisting parameter θ * given in (26) is chosen to minimize an upper bound on E θ T 2 γ th , we need to investigate whether θ * is close to the optimal unknown twisting parameter, that is the value that minimizes the actual value of E θ T 2 γ th . We consider the sum of two i.i.d Weibull RVs with shape and scale parameters equal to 0.5 and 1, respectively. In Fig. 8, we plot the upper bound (25) and the actual value of E θ T 2 γ th function in θ and for different value of the threshold γ th . We note that the exact computation of E T 2 γ th has a unique minimum which is closer to our choice θ * . Moreover, as γ th increases, the difference between the two minimizers becomes negligible and thus we tend to the optimal value. Another important deduction is that the second moment is slowly varying with respect to the twisting parameter θ especially in the neighborhood of the optimal value. Hence, our choice of θ * is actually reasonable since it almost results in approximately the biggest reduction of variance.

VI. CONCLUSION
In this paper, we developed an efficient hazard rate twisting technique for the estimation of the probability that a sum of independent RVs exceeds any given threshold. We presented a general procedure to find the best possible twisting parameter which leads to the possible largest reduction of the IS estimator variance for all possible values of the threshold. Besides, this approach, which seems to be consistent with the class of subexponential distributions, results in ensuring the asymptotic optimality criterion as the threshold goes to infinity. Numerical simulations showed that the optimized IS approach could reach the same accuracy as the naive MC simulation with a substantial computational gain. This alternative technique could serve as a benchmark to study the accuracy of future closed-form approximations of the quantity of interest.