On the Efficient Generation of α-κ-μ and α-η-μ White Samples with Applications

This paper is concerned with a simple and highly efficient random sequence generator for uncorrelated α-κ-μ and α-η-μ variates. The algorithm may yield an efficiency of almost 100%, and this high efficiency can be reached for all special cases such as α-μ, κ-μ, η-μ, Nakagami-m, Nakagami-q, Weibull, Hoyt, Rayleigh, Rice, Exponential, and the One-Sided Gaussian.This generator is implemented via the rejection technique and allows for arbitrary fading parameters. The goodness-of-fit is measured using the Kolmogorov-Smirnov and Anderson-Darling tests. The maximum likelihood parameter estimation for the κ-μ distribution is proposed and verified against true values of the parameters chosen in the generator. We also provide two important applications for the random sequence generator, the first one dealing with the performance assessment of a digital communication system over the α-κ-μ and α-η-μ fading channels and the second one dealing with the performance assessment of the spectrum sensing with energy detection over special cases of these channels.Theoretical and simulation results are compared, validating again the accuracy of the generators.


Introduction
In nearly all fields of science, simulation is a strikingly powerful tool widely adopted to help develop a better understanding of some phenomenon under investigation.Particularly in engineering, it is used, for instance, to successfully test equipment, algorithms, and techniques, and, to some extent and whenever applicable, to avoid or minimize timeconsuming, costly, and inexhaustible field trials.Wireless communications are no exception and in this challenging, lively, and unkind area, with systems becoming increasingly more complex, both industry and academy engage themselves in developing simulators.Such simulators for wireless communications almost certainly include a block for the fading channel.
The fading channel can be described by a number of models.Among them, the general models, namely, --, -- [1], and some particular cases such as - [2], - [2], and - [3], have been gaining wide acceptance .Their flexibility renders them adaptable to situations in which none of the traditional distributions yield good fit [2,3].In addition, their applicability has been recognized in practical and real scenarios.Field measurements carried out in diverse propagation environments have shown that, in many situations, these models better accommodate the statistical variations of the propagated signal [1, Section VII], [2,[7][8][9][10]26].In this sense, developing and ameliorating methods in order to simulate the -- and -- fading models and their special cases for arbitrary values of their parameters are of paramount importance.One first step in Section 6 the average error probability of the BPSK modulation over the -- and -- channels and the performance of the spectrum sensing over the - and - channels are presented.Some conclusions are drawn in Section 7.

Proposed Algorithm: Preliminary Results
In this section we present a preliminary proposed algorithm.However, in Section 4 the definitive and more efficient algorithm will be presented.
The majorizing hat function () used here is given as [30] where ℎ () is the majorizing density, , , and  0 are coefficients to be obtained for the specific fading model so that  () can majorize  --  () for all , and  -- P () is the desired probability density function (pdf) given in terms of the normalized envelope , with -- standing for -- or --.The parameter  is given in an exact form as where erf(⋅) is the error function.The coefficient  0 is obtained as the solution of (/)

𝑥-𝑦-𝑧 𝑃
() = 0.In all cases, the parameter  0 can be easily found numerically using wellknown software tools such as Mathematica or MATLAB.The coefficient  is found as the mode of the pdf; that is,  =  --  ( 0 ).Finally, the coefficient  is found as Algorithm 1 summarizes the steps for generating the desired sequences.The probability of acceptance in step 7 is 1/.U(0, 1) is the uniform distribution over the unit interval (0, 1].The rejection method is well known and it is described in detail in [31].Notice that the function ℎ() has the form of a truncated-Gaussian density.Random variables with pdf ℎ () can be generated in a fast and accurate way by truncated-Gaussian random variables generation methods (e.g., [32]).

2.1.
The -- Distribution.For a fading signal with envelope  and normalized envelope  = /  √E(  ), the normalized -- envelope pdf is given as [1] where  > 0 is a parameter describing the nonlinearity of the propagation medium,  > 0 is the ratio between the total power of the dominant components and the total power of the scattered waves,  > 0 is related to the number of multipath waves, and  ] (⋅) is the modified Bessel function of first kind and order ] [33].
(1) Define the distribution parameters (2) Find  0 , solving (/) In particular, the first derivative of the -- distribution is given by 2.2.The -- Distribution.For a fading signal with envelope  and normalized envelope  = /  √E(  ), the -- normalized envelope pdf is given as [1] where  ≥ 0 is the ratio between the in-phase scattered wave and the quadrature scattered wave and Γ(⋅) is the Euler Gamma function [33].
In particular, the first derivative of the -- distribution is given by

Numerical Results
In Figure 1, the empirical pdfs generated by the proposed method using 2 20 samples are contrasted with the theoretical density for different values of , , , and .In this figure, the solid lines correspond to the theoretical results whereas the symbols correspond to the generated random variates.The excellent agreement between theoretical and simulated results can be noticed.

Efficiency.
The acceptance proportion, or efficiency, is the performance measure of the acceptance-rejection method.It is the ratio between the number of samples accepted by the method and the total number of samples generated from the respective hat function (majorizing function).Figures 2  and 3 depict the efficiency curves for different values of the pdf parameters using the proposed algorithm.Hereafter, the solid lines correspond to the theoretical results whereas the symbols refer to the simulation results.The efficiency of the proposed method for the -- distribution is shown in Figure 2. Notice that the acceptance proportion decreases with the increasing of the parameter .As can be seen, the efficiency is rather small for  below 2 but increases rapidly as  reaches 2. The efficiency increases even further to reach almost 100% for  around 2.2.The well-known Nakagami- distribution is obtained by setting  → 0 and  = 2 in the -- distribution in which case  =  and the efficiency is around 75%, in agreement with [30].
The efficiency of the proposed method for the -- distribution is shown in Figure 3. Notice that the acceptance proportion increases with the increasing of the  parameter for a fixed  > 4. The Nakagami- distribution is obtained by setting  = 2, for  → 1, in which case  = /2, or, equivalently, for  → 0, in which case  =  in the -- distribution.In this case, once again, the efficiency stays around 75%, in agreement with [30].In all the cases, the acceptance ratio does not vary significantly with the variation of .In both, Figures 2 and 3, the acceptance proportion decreases with the increase of the fading parameter .This is a purely mathematical problem, in which we want to find a function (hat function) which is as close as possible to the distribution whose samples are to be generated.As it happens, the variations of the parameters provoke a change in the shape of the curves, which depart from that of the hat function, leading to an increase or a decrease in the efficiency.It is noteworthy that the samples are drawn with arbitrary parameters , , , and .Clearly, the acceptance proportion achieved using the proposed majorizing density is higher than when compared with a traditional uniform majorizing density.
A strikingly interesting result is shown next.Refer to Figure 2 for the efficiency of generating the -- random variable.It can be noticed that the efficiency reaches almost 100% for  ≅ 2.2.For instance, with  = 2.2,  = 3.5, and  = 4.25 the efficiency is 99.76%.A similar conclusion can be found for the efficiency of generating the -- random variable plotted in Figure 3.In this case the efficiency reaches almost 100% for  ≅ 3.5.
Figure 4 depicts the efficiency curves for the -- over different values of  and  with fixed  = 2.2.Notice that the efficiency starts above 85% ( = 0, - case with  = 2.2), increases, and decreases but is still above 95% for  > 1.
Figure 5 plots the efficiency curves for the -- over different values of  and  using  = 3.5.Notice that the efficiency starts above 95% ( = 1, - case with  = 3.5), increases, and decreases but is still above 92.5% for  > 1.

Goodness-of-Fit Test.
The difference between the theoretical and experimental distributions is minimal as visually perceived in Figure 1.However, in order to objectively quantify the performance of the random sequence generator for the -- and -- fading distributions, the Kolmogorov-Smirnov (KS) test is performed so that the empirical cdf and the hypothesized cdf are compared.As found in the literature (e.g., [29]), the measure of the fit accuracy is given by the -value.Table 1 reports -values obtained for the generated sequences with different values of , , , and .In all of the cases,  > 0.05, unveiling an excellent goodness-of-fit test result [31].It is well known that the Anderson-Darling test gives more weight to the tails than the KS test.Also, because the Anderson-Darling test is specific for the hypothesized distribution, this test is likely to be more powerful than the traditional KS test [31].For these reasons, we additionally have performed the Anderson-Darling goodness-of-fit test, with the objective of confirming the adherence of the generated random numbers also in the tails of their probability distributions.In general, critical values of the Anderson-Darling test statistic depend on the specific distribution being tested.We have tested two particular cases of the -- and -- distributions, since for these cases the exact critical values and, consequently, the  value are calculated analytically.The test has been performed for Weibull distribution (-- with  = 3.5,  = 0, and  = 1) and for the exponential distribution (-- with  = 1,  = 0, and  = 1).In the latter case the  value was 0.8352.In the first one the  value was 0.5236.These values reveal the adherence of the generated random numbers also in the tails of their probability distributions.

Main Result: A Near-100% Efficient Algorithm for Generating 𝛼-𝜅-𝜇, 𝛼-𝜂-𝜇 Variates, and Their Particular Cases
A modified procedure for a high-efficient and definitive algorithm can be noticed in this section.Let us consider first the -- density.From (12) of [3], it is possible to conclude that in the case of the -- distribution, for any set ( 1 , , ) and ( 2 , , ), the following holds: That is, given an -- distribution with parameters ( 1 , , ), another -- distribution with parameters ( 2 , , ), can be obtained by following the given transformation.In particular, knowing that an efficiency of almost 100% is achieved for -- samples with  = 2.2 (see Figure 4), a transformation of the kind   -- =  2.2 2.2-- can be used to attain any -- samples with this high efficiency.
Considering the -- distribution, for any set ( 1 , , ) and ( 2 , , ), the following holds: That is, given an -- distribution with parameters ( 1 , , ), another -- distribution with parameters ( 2 , , ) can be obtained by following the given transformation.Specifically, knowing that an efficiency of almost 100% is achieved for -- samples with  = 3.5 (see Figure 5), a transformation of the kind   -- =  3.5 3.5-- can be used to achieve any -- samples with this high efficiency.In such cases, the efficiency is kept constant and close to 100%, throughout the variation of the parameters.
Because - and - are particular cases of the -- and -- distributions (if we set  = 2), respectively, notice that for both, - and - distributions, the same high efficiency can be attained.In other words, in order to generate - samples with an almost 100% efficiency, the best choice is to generate -- samples with  = 2.2 and make the transformation  2  2-- =  2.2 2.2-- .For the - case,  2 2-- =  3.5  3.5-- .
In the same way, one can conclude that the high efficiency can be reached for all the particular cases of -- and -- distributions such as the well-known Nakagami-m, Rayleigh, and Weibull densities.All the particular cases of the -- and -- distributions can be found in [1,Section VI].
The steps for generating the desired sequences using the definitive algorithm are summarized in the Algorithm 2.

𝜅-𝜇 Random Variable Generation and Maximum-Likelihood Parameter Estimation
The generator ability in providing random samples following a given distribution can be alternatively verified by generating a large number of random variables and obtaining maximum likelihood (ML) estimates for the distribution parameters.In this section, as a particular case of the -- generator, we generated sample data sets of  independent identically distributed (i.i.d.) - random variables.Then, for each data set, we applied the ML parameter estimation and verified the estimates against true values of the generator parameters.

𝜅-𝜇
In particular, assuming  1 ,  2 , . . .,   as i.i.d.random variables and with [2, Equation (1)], the random samples have a joint pdf given by Equivalently, we can maximize the log-likelihood function (; ) = ln[ P (; )], which follows directly from (9) as to write From (11), one can see that it is necessary to simultaneously solve the following equations: in order to obtain κML and μML .
Taking the derivative of ( 10) with respect to , after some simplifications we have where  = 2√(1 + )  .In the same way Observe that  ] () in ( 14) depends on  with respect to both the order ] and the parameter .Hence, the derivative with respect to  in the last term of ( 14) is the sum of two terms, one only related to  and the other only related to ].As a result, we have Here  (1,0)   ] () is defined as the derivative of  ] () with respect to the order ] [33, Equation (9.6.42)]; that is International Journal of Antennas and Propagation where  0 (⋅) is the Digamma function [33, Equation (6.3.1)].
Fortunately, the function defined in ( 16) is available for direct usage in current numerical softwares, such as Mathematica, for instance, which makes ( 16) numerically tractable without additional difficulties.
In general, it is less computationally intensive to evaluate κML and μML iteratively by optimization algorithms, that is, finding ΘML that maximizes (; ) according to (11).
Here we use the iterative optimization algorithm  (⋅), available in the MATLAB software, to estimate the distribution parameters.In this case, we maximize the loglikelihood function applying the algorithm over the negative log-likelihood −(; ).However, we still make use of ( 13) and ( 14) in the estimator variance analysis.
The variance of an estimator is a measurement of its ability to perform reliably as it gives the degree of certainty in which the parameter is being estimated.In this context, the Cramér-Rao lower bound (CRLB) sets a lower limit for the variance of all unbiased estimators for  and gives the asymptotic variance for its ML estimator in a large sample size condition [34].Particularly, we can obtain the CRLB by evaluating the Fisher information matrix I  () contained in  random variables  1 ,  2 , . . .,   about the parameter  as [35] where the derivatives  ln   (; )/ and  ln   (; )/ are obtained, respectively, from ( 13) and ( 14) by setting  = 1.
One can readily verify that the derivatives are given by  ln   (; ) ln   (; ) The element [I()] 11 , for instance, is numerically solved as where   (; ) is the - pdf given by [2, Equation Also, it has been shown [35] that ΘML has an asymptotic multivariate Gaussian distribution with mean  and covariance matrix (1/)I() −1 .Thus, as  → ∞ Furthermore, it is straightforward to find [36] the confidence interval ((), ()) for ΘML = [κ ML , μML ], with confidence level of 95%, as We note that I() in ( 23) depends on , that is, the true value of the parameter, which is therefore unknown.However, ΘML converges in probability to  as  → ∞ [34], which makes it possible to infer that as  → ∞.As a result, for large sample sizes we can estimate I() by I( ΘML ) in ( 23) in order to compute the confidence interval.In this way we cover a useful range of the parameters  and , found in both indoor and outdoor multipath propagation environments [2,37].Taking advantage of notable properties of the ML estimation for large sample sizes, we calculate ML estimates of the parameters  and  for each sequence, according to (11), using the already cited algorithm  (⋅), available in the MATLAB software.The starting values of the estimates required by the algorithm are given as the true values of the parameters.

Performance of the 𝜅-𝜇
Figure 7 shows the sample mean of κML , (1/500) ∑ 500 =1 κML , against the true value of the parameter.In order to show the estimator variations about its sample mean, we also plotted in Figure 7 a confidence region defined by ±2 × (sample standard deviation of κML ), where the sample standard deviation is given by  [(1/500) ∑ 500 =1 κ2 ML − ((1/500) ∑ 500 =1 κML ) 2 ] 1/2 .We verify that, for a large sample size, κML is unbiased, from a practical point of view, for the useful range of .This is in accordance with the unbiased behavior that the ML estimators have in large sample size conditions [29].In addition, the confidence region of κML becomes broader as  increases; that is, the variance of κML increases with the value of the parameter, as observed from the results in Figure 6(a).
Similar results of the sample mean for μML against  are depicted in Figure 8.Likewise, μML is practically unbiased for the useful range of  and its variance increases with the value of the parameter, according to the results in Figure 6(b).
The unbiasedness of the ML estimators for large sample sizes, depicted in Figures 7 and 8, reveals the - generator ability to provide real random samples representative of a population with distribution   (; , ).This also alternatively confirms the excellent goodness-of-fit results given by the Kolmogorov-Smirnov and Anderson-Darling tests in Section 3, when applied to the -- generator.

Applications
In this section we give applications of the proposed random variable generators and use theoretical and simulation results for certifying the accuracy of these generators.

Average Error Probability of the BPSK Modulation over the 𝛼-𝜅-𝜇 and 𝛼-𝜂-𝜇 Fading Channels.
Here we analyze the bit error rate (BER) of the BPSK modulation over frequency-flat fading channels modeled by the -- and -- distributions.We assume coherent detection with matched filter or correlator receivers, for which the following vector channel model applies: the decision variable is  =  + V, where  is   the multiplicative fading with E[ 2 ] = 1,  = +1 represents a bit 1,  = −1 represents a bit 0, and V is the zero mean, additive white Gaussian noise (AWGN) with variance  0 /2.One possible analytical method employed for determining the performance of a mobile radio communication system is by evaluating the error probability as a function of a fixed signal-to-noise ratio (SNR) and then averaging the result over the probability density function of the SNR variations, which is governed by the particular envelope fading distribution.
International Journal of Antennas and Propagation For instance, the bit error probability of the BPSK modulation over the pure AWGN channel as a function of the received SNR  is given by where  =  2  b / 0 is the SNR for a particular value of the envelope ,  b / 0 is the ratio between the average energy per bit and the noise power spectral density, and erfc(⋅) is the complementary error function.Now, we must average  e () over the probability density function of ; that is, where
Applying a transformation of random variables, from (4) and ( 6), we have where  = E(  ) b / 0 is the average SNR when  is -- or -- distributed.
The integral in ( 26) was evaluated numerically.The average error probabilities curves for the -- and -- fading channels are presented in Figures 9 and 10.In these figures, the theoretical results (solid), from (26), and the simulated results (symbols), from a MATLAB program based on the vector channel model presented above, are plotted for the indicated fading conditions.It is clearly observed that, in spite of the variety of the fading conditions, excellent agreement between the estimated and theoretical results is shown, once more certifying the proposed algorithms for random variable generation.A myriad of different scenarios can be exercised for different values of the fading parameters.We omitted some particular cases (e.g., Rice or Nakagami-) for the sake of brevity.All the particular cases departing from the general models considered in this paper are in agreement with the particular cases presented in the literature (see, for example, the expressions presented in [38]).
It is worth mentioning that the application just described can be used to check the adherence of the generated random numbers to the tails of their probability distributions as follows.The agreement between theoretical and simulation results in the high  b / 0 regime is an evidence of a good adherence in the tail region, since this region governs the performance at high signal-to-noise levels.

Spectrum
Sensing over the - and - Fading Channels.Modern wireless communication systems are now facing a huge obstacle, spectrum scarcity.New services and applications appear every day, demanding increased bandwidth, new spectrum bands, or both.However, the currently adopted fixed spectrum allocation policy prevents those services and applications to be deployed in adequate pace.Nevertheless, recent studies have demonstrated that, in fact, the radiofrequency spectrum is quite underutilized in some areas and during some time [39].The cognitive radio (CR) [40] concept then came into scene, aiming at, among other things, opportunistic dynamic spectrum access to idle bands.In this situation, the network which owns the right of using the spectrum is called the primary network, and the cognitive radio network is usually referred to as the secondary network.
To detect the idle bands, also called spectral holes or whitespaces, the CRs must have some sort of spectrum sensing capability [41].Among the spectrum sensing techniques already developed, energy detection is one of the most attractive, since it has low implementation complexity and good detection power.In energy detection, a test statistic computed from the received signal energy and the noise variance is compared against a threshold so that the decision upon the occupation of the sensed channel is made.
Several studies consider the problem of spectrum sensing with energy detection over the pure AWGN channel, an approach, that is, by far unrealistic since typical wireless communication channels are also subjected to fading.Then, it is of paramount importance to access the performance of a spectrum sensing technique taking into account the channel fading.
There are several fading channel models available in the literature.Among them, two well-accepted models deserve attention due to their ability for accurately modelling several real channel conditions in practice.They are the - [2] and - [2] fading channel models.Both are special cases of the -- and -- fading models considered in this paper, if  = 2 is adopted.
In this section we apply -- and -- random variates, specialized to - and - variates, to analyze the performance of the energy detection over fading channels.

System Model.
The discrete-time model for the hypothesis test associated with the spectrum sensing problem is given by where H 0 denotes an idle channel state and H 1 denotes a busy channel,   is the th received signal sample collected by the CR,  = 1, . . .,  during the sensing interval,   is the zero mean Gaussian thermal noise sample generated at the receiver input,   is the primary transmitted signal sample, and ℎ represents the channel fading envelope, which is assumed to be constant during the sensing interval.From the received signal, the test statistic for the energy detector is computed according to where  2 =  0  is the thermal noise variance measured in the bandwidth , with  0 being the noise power spectral density.The number of samples  relates with the sensing time  and the bandwidth  through the time-bandwidth product  = , leading to  = 2.

Results
. The performance of a spectrum sensing technique is often measured in terms of the probability of detection,  d , and the probability of false alarm,  fa .When the primary network signal is present in the sensed channel,  d is the probability of declaring it indeed present.When the primary network signal is absent,  fa is the probability of declaring it present.A large value of  d translates into a small probability of interference from the secondary in the primary network.A small value of  fa translates into an increased throughput of the secondary network due to a more efficient use of spectral holes.These probabilities are often traded in a receiver operating characteristic (ROC) curve, which shows the values of  fa versus  d as the decision threshold is varied.
Figure 11 shows analytical (lines) and simulation (symbols) results of the energy detection over the - (-- with  = 2) and - (-- with  = 2) fading channels.We have considered  = 50 samples ( =  = 25),  = 3, 5, 8, and 10 dB.The simulation results were obtained from 100,000 Monte Carlo runs.The analytical results were obtained by numerically evaluating (6) and (11) of [42] in the case of -, and by evaluating (3) and (14) of [43] in the case of -.We have used the Mathematica software package to solve the above equations.The minimum and maximum values of the decision threshold (resp.,  min and  max ) are also reported to facilitate the reproduction of our results.Specifically, for  = 0.1 and  = 1,  min = 31, and  max = 100, for  = 1.5 and  = 0.8,  min = 30, and  max = 85, for  = 2 and  = 4,  min = 33, and  max = 86, and for  = 2 and  = 2,  min = 27, and  max = 80.
From Figure 11 the adherence between simulation and analytical results is apparent, validating, again, our random sequence generator.

Figure 11 :
Figure 11: ROC curves for - and - with different values of , , and .

Table 1 :
values of the KS test.