Discrete approximations of continuous probability distributions obtained by minimizing Cramér-von Mises-type distances

We consider the problem of approximating a continuous random variable, characterized by a cumulative distribution function (cdf) F(x), by means of k points, x1<x2<⋯<xk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_1<x_2<\dots < x_k$$\end{document}, with probabilities pi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_i$$\end{document}, i=1,⋯,k\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,\dots ,k$$\end{document}. For a given k, a criterion for determining the xi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_i$$\end{document} and pi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_i$$\end{document} of the approximating k-point discrete distribution can be the minimization of some distance to the original distribution. Here we consider the weighted Cramér-von Mises distance between the original cdf F(x) and the step-wise cdf F^(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{F}}(x)$$\end{document} of the approximating discrete distribution, characterized by a non-negative weighting function w(x). This problem has been already solved analytically when w(x) corresponds to the probability density function of the continuous random variable, w(x)=F′(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w(x)=F'(x)$$\end{document}, and when w(x) is a piece-wise constant function, through a numerical iterative procedure based on a homotopy continuation approach. In this paper, we propose and implement a solution to the problem for different choices of the weighting function w(x), highlighting how the results are affected by w(x) itself and by the number of approximating points k, in addition to F(x); although an analytic solution is not usually available, yet the problem can be numerically solved through an iterative method, which alternately updates the two sub-sets of k unknowns, the xi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_i$$\end{document}’s (or a transformation thereof) and the pi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_i$$\end{document}’s, till convergence. The main apparent advantage of these discrete approximations is their universality, since they can be applied to most continuous distributions, whether they possess or not the first moments. In order to shed some light on the proposed approaches, applications to several well-known continuous distributions (among them, the normal and the exponential) and to a practical problem where discretization is a useful tool are also illustrated.


Introduction
Using a discrete approximation of a continuous random variable (rv) is a procedure that is often adopted in many problems where uncertainty is present and needs to be taken into account.Substituting a continuous probability density function (pdf) with an approximating probability mass function (pmf), supported on a (possibly) finite number of points, can heavily reduce the computational burden required for determining a numerical solution for the problem at hand, and can produce an approximate solution whose degree of accuracy is still acceptable.This is particularly true when a problem involves several quantities that should be modelled as continuous rvs and the solution depends on some complex function thereof.The exact solution in this case could be obtained by applying some multivariate numerical integration technique whose computational cost may be severe and dramatically increasing with the number of rvs involved.Approximating each continuous rv through a properly chosen discrete rv allows the researcher to avoid numerical integration and resort to enumeration, which is much easier to manage (Luceno 1999).
An interesting application of discrete approximation of continuous distributions can be found in the field of insurance.If one is required to determine the distribution of the total claim amount S = N i=1 X i corresponding to a random number N of i.i.d.claims whose size X i is assumed to follow some continuous probability distribution, rather than resorting to integral convolution to determine the exact distribution of S, one can approximate the distribution of the claim size X i through discretization and then apply Panjer's formula (Panjer 1981), which is a recursive formula that exactly calculates the distribution of the total, holding when the claim size is arithmetic with span size h > 0 and the distribution of the number of claims N belongs to the (a, b, 0) class.
In the field of quantitative finance, a similar application is the following.Let L L L = (L 1 , . . ., L d ) denote a vector of possibly dependent rvs, each one representing a loss on a particular trading desk, portfolio or operating unit within a firm, over a fixed time period.Sometimes we need to aggregate these losses into a single rv, typically the sum L + = d i=1 L i , on which one can calculate a measure of the aggregate risk, for example, the Value-at-Risk, which is nothing else than the quantile at a prespecified level 0 < α < 1 of the distribution of L + .However, determining this measure of risk, when the joint distribution of L L L is fully specified, requires the computation of the distribution of L + , which is not straightforward to derive even if all the rvs are mutually independent.A possible answer is represented by the discretization of the L i and the construction of a joint pmf approximating the joint distribution of L L L, on which the evaluation of the VaR of the aggregating function is much more straightforward (see, for example, Jamshidian and Zhu (1996), where the authors consider the multivariate Gaussian distribution).
In reliability engineering, the stress-strength model describes a system with random strength which is subject to a random stress during its functioning, so that the system works only when the strength is greater than the stress.The probability that a system correctly works is termed reliability (Johnson 1988).Evaluating the reliability of a system thus requires the knowledge of the probability distributions of both stress and strength; when these latter depend on several stochastic factors, the probability distributions (and then the reliability) are often not analytically tractable and hence some form of approximation is needed.Given a known functional relationship between stress (or strength) and its random subfactors and assuming the subfactors of the stress (or strength) are independent, one feasible approach of approximating the probability distribution of stress (or strength) is through discretization of the subfactors (see, e.g., English et al. 1996).
For an assigned continuous probability distribution and a fixed number of approximating points k, the matter is how to build an "optimal" discrete approximation.Several criteria have been used so far; here is a rough classification into four main categories: -moment equalization or moment matching: the discrete approximation is the one preserving as many moments as possible of the original distribution.Moments' matching is carried out through a procedure known as Gaussian quadrature (Golub and Welsch 1969); in its more authentic form, the support points of the discrete approximation and their probabilities are derived simultaneously.This is by far the oldest and most popular discretization technique; as stated by Miller and Rice (1983), "Few people would accept an approximation that did not have roughly the same mean, variance, and skew as the original distribution".However, its application is limited by the finiteness and existence of a closed-form expression for the first integer moments of the continuous rv: many random distributions commonly used in quantitative finance, for example, do not possess even lowerorder moments.Several variants of Gaussian quadrature have been proposed; for example, Tanaka and Toda (2013) suggested that the k support points of the approximating distribution have to be chosen "a priori" and the probabilities have to be derived by maximizing the relative entropy to an assigned "reference distribution", under the constraint that it matches as many moments as possible.Convergence properties of this maximum entropy method and applications to stochastic processes were later presented in Tanaka and Toda (2015), Farmer and Toda (2017); -preservation of the distribution function: the discrete approximation preserves, at each support point, the value of the cdf or, alternatively, of the survival function (Roy and Dasgupta 2001).Actually, this technique has been employed for constructing discrete counterparts of continuous distributions, by defining its pmf as p(i) = F(i)− F(i −1), with i integer; such construction automatically supplies a valid pmf, which preserves the cdf of the continuous distribution at any integer value i (Chakraborty 2015).Straightforward modifications have been proposed in order to handle finite supports consisting of possibly non-integer points; -minimizing the mean squared error between the assigned continuous rv X ∼ F and its approximation X , i.e., minimizing E F (X − X ) 2 , where the expected value is computed with respect to the distribution of X ; this yields the so-called optimal quantization, which is a well-known technique in signal theory (Gray and Neuhoff 1998).The minimization problem can be solved iteratively, by alternately solving two sets of conditions, one expressing each support point (also named "quantum") as the center of mass of the intervals of a partition of the support of the continuous rv, and the other expressing the endpoint of each of these intervals as the midpoint of the segments between successive quanta (Lloyd 1982); -minimizing a distance between the two distribution functions: the discrete approximation is obtained as the distribution minimizing some statistical distance between the two (continuous and step-wise) cdf.In Kennan (2006), an analytical solution is found for the optimal k-point discrete approximation when employing a particular class of distances.
In this work, we will concentrate on the last class of discrete approximations.More precisely, we propose the use of three distances: the Cramér-von Mises, the Anderson-Darling, and the Cramér distance between the cdf of the assigned continuous rv and the step-wise cdf of its k-point discrete approximation.They can be all regarded as weighted Cramér-von Mises distances with a proper choice of the weighting function.We will derive for these three cases the optimal solution (i.e., the k-point discrete approximation leading to the minimum distance) either analytically or computationally, providing in the latter case some details about the numerical procedure to be implemented.To the best of our knowledge, these distances have not been used regularly in the existing literature for the discrete approximation of a continuous rv.The aim of this paper is to shed some light on their use and explain their pros and cons.
The rest of the paper is structured as follows.In the next section, we will discuss distances between cdf in general and then present and solve the main research question, that is, finding an optimal k-point discrete approximation to a continuous random distribution through minimization of a statistical distance.The three statistical distances mentioned above will be considered and the corresponding solutions to the problem will be described and compared, also practically referring to some known parametric families of continuous distributions.Section 3 illustrates two applications of discretization to the (approximate) determination of the distribution or of some parameter of an assigned function of independent rvs, which is carried out by using each of the discrete approximations previously described and also other existing techniques.A software implementation in the R programming environment is presented in Sect. 4. Some comments and remarks are provided in the last section.mum (or, better, supremum) absolute difference between the two cdfs.The KS distance is commonly employed when one has to test whether an i.i.d.sample (x 1 , x 2 , . . ., x n ) is consistent with some known cdf F; in this case, the distance to be computed is between the assigned F and the empirical cdf F(x), defined as Another wide family of distances between cdfs is the following where w(x) is some non-negative function on R.
If we set w(x) ≡ 1, then we obtain the second-order squared Cramér distance (Cramér 1928), which is closely related to the so-called energy distance (Rizzo and Székely 2016).If we set w(x) ≡ f (x), then we obtain the so-called Cramér-von Mises distance, which can be also conveniently written as then we obtain the Anderson-Darling distance; with respect to Cramér-von Mises, this distance puts more weight in the tails of the distribution, i.e., where F(x) is close to either zero or one.Hanebeck and Klumpp (2008), also referring to the less explored multivariate case, remark "how statistical distances are used for both analysis and synthesis purposes.Analysis is concerned with assessing whether a given sample stems from a given continuous distribution.Synthesis is concerned with both density estimation, i.e., calculating a suitable continuous approximation of a given sample, and density discretization, i.e., approximation of a given continuous random vector by a discrete one".In this work, we focus on the latter aspect of research.
If one is interested in building an optimal k-point discrete approximation of a given continuous cdf F, i.e., a discrete probability distribution somehow resembling the original one, then he/she can find it as the discrete distribution minimizing, over all the k-point discrete distributions, the distance (1), computed between the cdf F and the cdf F of the discrete approximation, for a given choice of the weighting function w(x).
The problem can be more formally stated as follows.Suppose we want to approximate the continuous probability distribution F by a discrete probability distribution k and k i=1 p i = 1).Let us gather the support points and the probabilities of the discrete distribution in a vector η η η = (x 1 , . . ., x k , p 1 , . . ., p k ).Then, the optimal discrete distribution can be defined as the one (univocally identified by η η η) minimizing d(F(x), F(x; η η η)): η η η = arg min d(F, F; η η η).For the case w(x) ≡ f (x), corresponding to the Cramér-von Mises distance, an analytical solution is available (Kennan 2006).When w(x) is a piece-wise constant function, the problem was solved numerically by Schremp et al. (2006), through an iterative procedure based on a homotopy continuation approach.For any other possible choice of the weighting function w(x), we can solve the minimization problem (at least) numerically.
We start our analysis from the Cramér-von Mises distance.

Cramér-von Mises
The Cramér-von Mises distance between continuous cdf is one of the distinguished measures of deviation between distributions (Cramér 1928;von Mises 1931).For a probabilistic interpretation, see Baringhaus and Henze (2017).It is obtained by setting w(x) ≡ f (x) in (1), so that the distance between F and F becomes where F −1 denotes the mathematical inverse of the cdf F of X , which we assume to be strictly increasing over the support of X ; the distance in (2) is a particular case of The best k-point discrete approximation (that is, the one yielding the minimum value of the distance d r ) has been proved (Kennan 2006) to have, for any r > 0, equallyweighted support points x i (i.e., with probability 1/k each) such that F(x i ) = 2i−1 2k or, equivalently, Here we report the proof, fixing and adding some points left out in Kennan (2006).
Let us introduce the quantities q i and Q i , defined in the following way: q i = F(x i ), i = 1, . . ., k, setting q 0 = 0 and q k+1 = 1; Q i = F(x i ), i = 1, . . ., k, setting Q 0 = 0; it follows that Q k = 1.Therefore, the q i represent the values of the cdf of the continuous rv at the support points x i of its discrete approximation; the Q i represent the values of the cdf of the discrete approximating rv at its support points 1).
The distance (3) can be thus rewritten as and the first-order condition on the Q i gives from which Q i = (q i + q i+1 )/2; the first-order condition on the q i provides Fig. 1 Meaning of q i and Q i when constructing a k-point discrete approximation based on the minimization of the Cramér-von Mises distance (here k = 4).Since q i = F(x i ) and Q i = F(x i ), it is possible to alternatively represent any k-point discrete approximation as a set of k points (q i , Q i ) belonging to the the unit square from which Combining these two last results together, we obtain By substituting this expression in that for the q i , we obtain and then Note the particular form of the optimal solution: the specific continuous distribution enters the equation of the support points through the inverse of its cdf F, which is applied to the q i , which are "distribution-free" quantities; the probabilities p i are independent from F as well, being all constant.
Another interesting property of this solution can be pointed out.Letting F k be the set of discrete distributions with k support points, we have just proved that for each r > 1 the k-point discrete approximation F ∈ F k that satisfies inf is the discrete uniform distribution on the set of point recalling the relationship between the L r norm and the supremum norm (see e.g.Stein and Shakarchi 2011), one obtains and deduces from this that inf that is, for a fixed integer k, the best approximating distribution obtained by minimizing d r is the same we would obtain by minimizing d K S .
Another peculiarity of the optimal solution or, better, of the statistical distance used as a criterion for finding an optimal solution, is its all-encompassing applicability, since it does not require the existence and finiteness of any integer moment of the original continuous distribution: the distance (3), since it can be rewritten as an integral over the unit interval of a quantity which is finite, can be always computed and always possesses a global minimum, corresponding to the solution derived above.Therefore, it is possible to derive such a discrete approximation for Student's t, say, with any value of the degree-of-freedom parameter, and other heavy-tailed distributions, which would not be possible in general if using the moment-matching or quantization techniques (comprising the method by Drezner and Zerom (2016), which can be considered as a compromise between the two techniques): the finiteness of the first 2k − 1 moments is required by the former, the finiteness of the first two moments by the latter.It can be also shown (Barbiero and Hitaj 2022) that lim k→∞ Fk (x) = F(x) ∀x ∈ R, i.e., the approximating k-point discrete rv converges in distribution to the original rv.
The distance (3) can be further rewritten as the sum of k + 1 integrals, and its minimum value is thus equal to since for the optimal solution q Furthermore, as one can expect, the minimum distance is a decreasing function of k for any fixed r > 0: by increasing the number of points, the optimal approximating discrete distribution gets closer (in terms of Cramér-von Mises distance) to the continuous one.Table 1 displays, just for illustrative purposes, k-point approximations (k = 5; 6; 7) of a standard normal and an exponential distribution (with unit rate parameter).For each k and for both distributions, values of expectation and variance (simply indicated as μ and σ 2 ) are reported, in order to be easily compared with the analogous values of the continuous distribution.
For the standard normal distribution, variance and kurtosis of the k-point discrete approximation monotonically converge to the corresponding value of the original continuous distribution rather slowly: kurtosis, in particular, when k = 100, is 2.834.This result could have been expected: although the cdf of this discrete approximation converges to the cdf of the original continuous rv, this does not mean that for a finite k the two functions are so similar: for the normal case, recall that the pdf of the continuous rv is displayed through the classical bell-shaped curve, whereas the pmf of the discrete approximation has uniform probabilities, and this overall translates into a mismatch of (even) moments.Recall that moment equalization would produce discrete distributions with a very large range and very small probabilities assigned to the extreme support points (Barbiero and Hitaj 2022).
For the discrete approximation of the exponential distribution, expected value, variance, skewness, and kurtosis tend monotonically and asymptotically to the value of the parent distribution (1, 1, 2, and 9, respectively), but for finite k, it underestimates kurtosis to a large extent; when k = 7, the value of kurtosis for the discretized exponential distribution is 2.779 (its expected value is 0.951, its variance 0.686 and its skewness 0.963); when k = 100, it is 6.662 (its expected value is 0.997, its variance 0.960 and its skewness 1.759).

Anderson-Darling
What happens if we consider a weighted Cramér-von Mises distance, for example the Anderson-Darling distance, characterized by the weighting function w We may expect that the best approximating distribution has now unequal probabilities or that the support points are no longer quantiles of   Minimizing the Anderson-Darling distance is equivalent to minimizing the following quantity, which consists of the sum of k + 1 contributions: with respect to the q i and Q i .The first-order condition for Q i is The first-order condition for q i implies from which we obtain again The solution derived by Eqs.( 6) and ( 7) cannot be expressed in an analytic closed form for each q i and Q i .However, a simple iterative algorithm can be implemented with the aim of recovering their values numerically; as initial guess values for q i and Q i , we adopt their optimal values found by minimizing the Cramér-von Mises distance.In a similar fashion as done in Pavlikov and Uryasev (2018), one can think of alternatively updating the values of the probabilities p i and the values of the discrete points x i (or, better, their probability transforms q i ) till convergence.Note that in Pavlikov and Uryasev (2018) all the p i are initially set equal to 1/k, as for the analytical solution based on the Cramér-von Mises distance.The algorithm works as follows: or any large positive value) and max = 10 −6 (or any arbitrarily small positive value, to be used for checking convergence of the solution) 3.While (t) > max : (a) Update the iteration index t ← t + 1 (b) Update the q i according to (7): (d) Derive the updated probabilities for the discrete rv: i−1 (e) Calculate the maximum absolute deviation between two consecutive iterations in terms of p i : Alternatively, other distances can be used to compare the probability vectors obtained in two consecutive iterations, as the Euclidean distance The solution is expressed in terms of q i (probability transformation of i and p

(t) i
The number of iterations required clearly depends on the threshold max and on the number of points k: diminishing max or increasing k leads to a larger number of iterations; for plausible values of (say 10 −8 ) and k (say smaller than one hundred), the computation times are in the order of fractions of a second.By numerical inspection, the Anderson-Darling weighting function leads to an inverted U-shaped trend for the p i , which turn out to be symmetrical around the central value or values, according to whether k is odd or even, i.e., p j = p k− j+1 , j = 1, . . ., k; see Table 2.As for the optimal q i values, they present a form of symmetry around the central value(s), such that a continuous symmetrical distribution remains symmetrical after discretization (see Tables 2 and 3a, referring to the discretization of a standard normal rv).
If we consider the exponential distribution with unit rate parameter and we derive the k-point discrete approximation by minimizing the Anderson-Darling distance with k = 7 (see Table 3b), its expected value turns to be equal to 0.961, its variance 0.743, its skewness 1.147, and its kurtosis 3.424.Although not so close to the values of the underlying continuous distribution, the discrepancies are smaller if compared to those resulting from Cramér-von Mises approximation (see the previous subsection).
We note that using a distance d(F, , with any other possible positive weighting function w(x), even asymmetrical, supplies (although in general only numerically) values of q i and Q i that do not depend on the specific cdf F(x); the original distance can in fact be rewritten as k i=0 q i+1 q i (t − Q i ) 2 w(t)dt and the first-order condition on the Q i becomes q i+1 q i (t − Q i )w(t)dt = 0 and thus the Q i that solve this condition can be expressed as a function, dependent on w(x), of the q i and q i+1 only.The simplest and most apparent case is provided by the "unweighted" Cramér-von Mises distance discussed in Sect.2.2, where such values are available through simple analytical formulas involving just the index i and the number of support points k.

Cramér distance
We now consider the Cramér distance between the cdf of a continuous rv and that of its discrete approximation: which can be also rewritten as This last expression highlights how, differently from the two distances previously examined, the Cramér distance may not be always computed: for some cdf F, the first and/or the last term of the sum, in fact, may be not finite: although the two corresponding integrals are computed over a limited interval, the integrand function may be not finite (when t tends to 0 or t tends to 1, the denominator F (F −1 (t)) may tend to zero; and in the former case, for example, it may tend to zero as fast as t 3 or faster and then the improper integral would not converge).We will see an example later.
Assuming that the Cramér distance can be computed, the first-order condition on the q i leads again to while the first-order condition on the Q i leads to for i = 1, . . ., k.This condition therefore depends on the specific distribution F of the continuous rv X .
We will now analytically develop condition (9) for several choices of F, by examining some well known parametric families.We will see that in general Eq. ( 9), combined with the condition on the q i , does not lead to a closed-form solution of the discrete approximation, which must be solved numerically, for instance resorting to the iterative procedure sketched out in Sect.2.3.

Normal
In case of a standard normal rv, with cdf e −t 2 /2 dt, for which the following equality holds (see, for example, Owen (1980), formula 1000): where φ(x) = Φ (x), the first-order condition on the Q i becomes which, along with the first-order condition on the q i , allows us to determine the optimal solution numerically, following analogous steps to those sketched out in Sect.2.3.If instead of considering a standard normal rv, we focus on a generic normal rv with parameters μ and σ 2 , with cdf F(x) = Φ( x−μ σ ), the first-order condition on the Q i would not change; both the left and right members of (9), in fact, are simply multiplied by σ .Table 4 displays the values (x i , p i ) of the best k-point approximation of a standard normal rv, for k = 5; 6; 7. We note that, as for the discrete approximation based on the minimization of the Anderson-Darling distance, the support points are symmetrical around zero and the probabilities satisfy p j = p k− j+1 , j = 1, . . ., k.

Exponential
In case of an exponential rv, with cdf F(x) = 1 − e −λx , x > 0, and quantile function we have that the first order condition on the Q i can be written as Table 5 displays the values (x i , p i ) of the best k-point approximation of an exponential rv with unit parameter for some values of k.It is very important to notice that empirical inspection shows that this type of approximation is characterized by a decreasing pmf, which thus resembles the decreasing trend of the pdf; on the contrary, the discrete approximation derived from the minimization of the Anderson-Darling distance possesses an increasing-decreasing pmf; the discrete approximation based on the minimization of Cramér-von Mises distance possesses a constant pmf.If we

Cauchy
In case of a standard Cauchy rv, with cdf F(x) = 1 π arctan x + 1 2 , x ∈ R, we have that the first-order condition on the We highlight that it is possible to find a k-point discrete approximation for the Cauchy distribution by minimizing the Cramér distance (as well as Cramér-von Mises or Anderson-Darling), although such a distribution does not possess any positive integer moment and then any discretization technique based on some form of momentequalization -among others, moment matching, which requires equalization of the first 2k − 1 moments, but even the technique presented by Drezner and Zerom (2016) -would be not applicable at all.

Logistic
If we consider the standard logistic distribution with cdf , and the first-order condition on the Q i can be written as: − log q i 1 − q i and then We note that the latter condition is the same as condition (6), obtained for the optimal k-point approximation based on the minimization of Anderson-Darling distance.This occurs since for the logistic distribution the equality f (x)[F(x)(1 − F(x))] −1 = 1 holds for any x ∈ R, and hence the Cramér distance coincides with the Anderson-Darling distance.
Similarly to what happens with the normal distribution, if we consider a nonstandard logistic distribution with location parameter μ and scale parameter σ , with cdf F(x; μ, σ ) = and quantile function F −1 (u; μ, σ ) = μ + σ log u 1−u , it is straightforward to see that the first-order condition on Q i would remain the same.

Lomax or Pareto
For the Lomax distribution with scale parameter λ > 0 and shape parameter α > 0, the expression of the cdf is 123 from which In particular, if α = 2, formula (12) reduces to the following expression If α = 1, the first-order condition on the Q i can be written as , from which If α ≤ 1/2, it can be proved that the Cramér distance cannot be computed; in fact, for any feasible value of λ and k, since the integrand function of the last integral in (8), )dt, turns out to be proportional to (1 − t) 1−1/α , then, being α ≤ 1/2, the (improper) integral is not finite for any feasible value of q k .We highlight that the Lomax distribution does not possess any integer moment for α ≤ 1, so the proposed procedure is still able to produce a k-point discrete approximation even if the distribution does not possess a finite expectation, for 1/2 < α ≤ 1.
We would obtain the same results if we considered the Pareto rv with cdf F(x) = 1 − λ x α , x > λ, by simply substituting x with λ + x.
In general, it is not possible to obtain a closed-form optimal solution, unless for small values of k.For example, if λ = 1 and α = 2, for k = 2, the conditions on the q i are q 1 = Q 1 /2 and q 2 = (1+Q 1 )/2; the condition on Substituting the first two expressions for q 1 and q 2 into the last condition, we obtain, after a few simple passages, that 1 4 and then a second-order equation in Q 1 whose unique feasible solution is Q 1 = 2 3 .Consequently, q 1 = 1 3 and q 2 = 5 6 ; the best discrete approximation consists of the points x 1 = 3 2 − 1 and x 2 = √ 6 − 1 with probabilities p 1 = 2/3 and p 2 = 1/3.

Power function
The cdf of a power function rv is The first-order condition on Q i turns into If c = 1 (corresponding to the case of a uniform distribution) then Q i = (q i +q i+1 )/2 and we obviously obtain the same solution as in Sect.2.2.If c = 2, then we have thus the cumulative probability Q i of the optimal solution corresponds to the arithmetic mean of the two consecutive probabilities q i and q i+1 and their geometric mean.
It can be numerically shown that for values of c larger than 1 (when the pdf is increasing), the k-point discrete approximation has an increasing pmf; on the contrary, for 0 < c < 1 (when the pdf is decreasing), the k-point discrete approximation has a decreasing pmf.
In general, it is not possible to derive the optimal solution in a closed-form, unless for small values of k.For example, if b = 1 and c = 2, for k = 2, the conditions on the q i are q 1 = Q 1 /2 and q 2 = (1+Q 1 )/2; the condition on Q 1 is Q 1 = 1 3 (q 1 +q 2 + √ q 1 q 2 ).Substituting the first two expressions for q 1 and q 2 into the last equation, we obtain, after a few simple passages, that 4Q

The case k = 1
We have always implicitly assumed so far that the number of points k by which one approximates the assigned continuous random distribution is greater than 1.It is immediate to realize that if one wants to approximate the distribution of the rv X through one value only, then, except for the cases where the selected statistical distance cannot be computed (see Sect. 2.4.5) the "optimal" approximating value x 1 is the median of F, F −1 (0.5), since the only effective condition to be satisfied (for all the three distances examined) would be q 1 = (Q 0 + Q 1 )/2 = 1/2, being Q 0 = 0 and Q 1 = 1.We note that optimal quantization (Lloyd 1982) would return the expected value E F (X ), as far as it exists and it is finite, as the optimal approximating value of F.

Location-scale transformation
Let us consider a continuous rv X ∼ F X and its optimal k-point approximation derived from the minimization of the Cramér-von Mises or Anderson-Darling distance, We know that the optimal q * i and Q * i are only functions of k and of the selected distance, but do not depend on the specific F X .Then consider the location-scale transformation Y = a +bX, a ∈ R, b > 0. Since for such a transformation we have F Y (y) = F X ( y−a b ) for any real y and F −1 Y (u) = a +bF −1 (u), 0 < u < 1, the optimal k-point approximation of Y (derived by using the same distance as for X ) is represented by = a + bx i and the same p * i as before, which means that the location-scale transformation applies also to the discretized rv.
This property does not hold in general for the discrete approximation based on the minimization of the Cramér distance: the optimal values q * i and Q * i now depend on the specific form of the cdf, and F Y may not belong to the same family as F X , unless the cdf F X is a location-scale family of distributions itself: in this case, the property still holds, as the condition ( 9) does not change if we consider a location-scale transformation of the cdf (recall the examples with the normal and logistic distributions in Sect.2.4).

Example of application
Often a researcher in the statistical field is required to determine the distribution (or parameters) of some complex function of several (independent and continuous) rvs T = t(X 1 , X 2 , . . ., X d ).Multidimensional integration techniques should be employed, but often, due to the complexity of t and to the high dimensionality d, they are either cumbersome to apply or not applicable at all.One can then resort to Monte Carlo simulation, i.e., simulating (independently) a huge number N of pseudorandom values from X 1 , X 2 , . . ., X d and then calculating the corresponding N values of the transformation t, which can be regarded as a random sample drawn from T .An alternative is to find an approximate solution via approximation-by-discretization and enumeration.This consists of substituting each X i with a discrete approximation Xi and then determine the corresponding pmf of T = t( X1 , X2 , . . ., Xd ) by "enumeration", based on the joint pmf of ( X1 , X2 , . . ., Xd ) over the Cartesian product of the single supports: In order to illustrate and compare the discretization techniques proposed in Sect.2, in the following subsection we will consider a case where it is possible to derive the exact distribution of T analytically, which allows us to evaluate the statistical performance (in terms of degree of approximation, to be measured through some index) of each technique; in Sect.3.2, we will analyzize a more complicated case where the exact distribution of T is not analytically computable, but a parameter of interest can be recovered numerically.Approximation is obtained by applying one of the univariate discrete approximations of Sect. 2 to each rv independently and then reconstructing the distribution of the sum by "enumeration"

Sum of exponential random variables
).However, let us approximate each X i through one of the discrete approximations we illustrated, and then reconstruct the approximated distribution of the sum, FS .We can then evaluate the degree of approximation by using a measure of discrepancy between F S and FS ; we can use the KS distance, rather than one of the statistical distances employed for the univariate approximation.Here we consider the sum of d = 3 exponential rvs; λ i = 1/2, α i = 1, i = 1, 2, 3, for which we know that S ∼ Gamma(3, 1/2).Table ( 6), for different values of k (from 5 to 11), reports the KS distance between the exact and the approximated cdf of S, obtained according to the three univariate discrete approximations for the X i .It is easy to see how the KS distance decreases with k for all the methods: by increasing the number of approximating points for the X i , it is legitimate to expect that the discrepancy between the true and the approximate cdf decreases, independently from the measure used.Moreover, we notice that the discrete approximation based on Cramér distance overperforms the other two approximations, and its relative level of accuracy quickly improves with k: for k = 11, the KS distance provided by the approximation based on the minimization of Cramér distance is much less than one half of the KS distance provided by the approximation based on the minimization of the Cramér-von Mises distance.As we already remarked, the discrete approximation of the exponential distribution based on the minimization of the Cramér distance seems to resemble the shape of the continuous pdf better than the discrete approximations based on the Cramér-von Mises and Anderson-Darling distances, and this also positively reflects when calculating the approximate distribution of the sum of three i.i.d.exponential rvs.
The graph in Fig. 2 displays the curve of the true cdf of the sum, and the three stepwise cdf of the three discrete approximations.It is quite apparent that the latter struggle in resembling the continuous curve in its right tail: it is visible to an unaided eye, in fact, for relatively high values of x (say, between 10 and 15) the three approximations show their maximum absolute error (which corresponds to the value of KS distance).
, where M is the applied torque, t is the wall thickess, W is the width, and H the height of the tube.Let X , M, t, W , and H be mutually independent rvs.Consider the following parametric set-up: 3, σ H = 0.03) (with N denoting the normal distribution) which accord with that in Roy and Dasgupta (2001).Let the standard deviation of the normal strength X be 60, and assume that the mean of X can vary from 520 to 970 units by steps of 10, thus generating an array of n = 46 scenarios.
Although the exact evaluation of R is unfeasible here, as suggested by a referee, its value can be obtained numerically as follows.Denoting by φ μ,σ and Φ μ,σ the pdf and cdf, respectively, of a normal rv with mean μ and standard deviation σ , putting for (t, m, w, h) ∈ R 4 , the reliability parameter R can be expressed as and can be approximated by where Q is the hyper-rectangle defined as and γ t , γ M , γ W , γ H are positive scale factors to be chosen suitably large.The computation of R Q is easily done by using the function cuhre that is included in the package cubature (Narasimhan et al. 2022) of the statistical software environment R (R Core Team 2022).
Alternatively, one can resort either to Monte Carlo simulation or to the approximation-by-discretization approach, by using one of the methods described in the previous sections.
We will regard the value of reliability recovered by numerical evaluation through the cuhre function (with all the γ scaling factors set equal to 9) as the actual (true) value and thus indicate it simply with R; the value obtained through Monte Carlo simulation is denoted by RMC ; the values obtained by the approximation-by-discretization approaches are denoted by R, with a superscript identifying the specific discretization: GQ (Gaussian Quadrature), DZ (Drezner and Zerom 2016), CvM (Cramér-von Mises), C (Cramér), AD (Anderson-Darling).In order to measure the degree of accuracy of each technique, one can consider the following synthetic measures, which are: Mean Deviation (MD), Mean Absolute Deviation (MAD), Root Mean Squared Deviation (RMSD), defined as follows, where the subscript i refers to the i-th scenario.
We considered a number N = 10 millions of pseudo-random simulations for the Monte Carlo approach, and a number k = 5 of approximating points for the discretization approaches.
The results, reported in Table 7, show that Monte Carlo simulation provides overall the best results in terms of MAD and RMSD; the proposed discretization methods, based on the minimization of the Anderson-Darling and Cramér distances, perform better if compared to the Gaussian quadrature method in terms of MD, MAD, and RMDE; they perform worse than the method proposed in Drezner and Zerom (2016), which however requires a much considerable computational effort required by its inner numerical minimization routine and has a narrower applicability, as already remarked in Sect.2.2.We remark that in this example we focused just on discretization techniques which are somehow comparable: according to their specific criterion, they all compute both the k support points and the corresponding k probabilities simultaneously.Other techniques, mentioned in the Introduction, such as the maximum entropy method, first assign the support points a priori and then calculate the probabilities only.

Software implementation
Code in the R programming environment has been developed, which implements the different routines used for finding the optimal discrete approximations.For the Cramér distance, in particular, several parametric distributions are considered.Some functions are also available for plotting discretized distributions, which makes graphical comparison to the original continuous distribution more effective.
For the Cramér-von Mises and Anderson-Darling distances (and potentially, for any other distribution-free distance), the function Discr has to be used, whose arguments Note that the values of the mean of X , μ X , have been actually chosen in order to induce a range of values for R roughly spanning the entire interval between 0.01 to 0.99 are the number of points k and the type of distance (CvM and AD).It returns a list containing the vectors of p i and q i (along with the vector of Q i ), from which one can extrapolate the support values x i through the quantile transformation of the assigned distribution function F.
For the Cramér distance, the main function is called DiscrF; its arguments are the number of points k by which we want to approximate the continuous random distribution; the type of continuous distribution (a string identifying it), along with the value of its parameters (a vector, par).At the moment, a few families of continuous distributions can be selected, namely, those discussed in Sect.2.4, for which the first-order condition on the Q i is available in a closed-form: normal (norm), exponential (exp), Cauchy (cauchy), logistic (logis), Lomax (lomax), power function (power).The output of the function is a list, containing the vector of probabilities p i and the vector of support values x i of the discrete approximation, along with the vectors of q i and Q i .
The companion function moments receives as input a vector of support points and a vector of corresponding probabilities possibly obtained as a result from DiscrF or Discr, and computes the expectation, variance, skewness and kurtosis of the related discrete rv.This is useful if one wants to keep under control the effects of discretization over the first (normalized) moments, since we know that the techniques we introduced do not offer any guarantee, in general, that the moments of the continuous rv are preserved.
The graphical function plotdist receives as its first argument the result of Discr or DiscrF and plots the corresponding discrete approximation (its pmf or its cdf, to be selected through the argument plot, to be set equal to "pmf" or "cdf").This can be plotted over the unit square (by setting the argument xaxis equal to "q": on the x axis the q i , on the y axis the p i or the Q i , see Fig. 1); or on the usual R × [0, 1] space (by setting the argument xaxis equal to "x": on the x axis the x i , which should be supplied if the first argument comes from Discr; on the y axis the p i or the Q i ).
In the supplementary material, available at https://tinyurl.com/STPA-D-21-00382, the relevant R code along with some examples is supplied.

Conclusion
In this work, we discussed a class of discretization techniques that calculate a k-point discrete approximation of a continuous rv by minimizing a distance between the two cdf.For one distance in this class, the optimal discrete approximation (i.e., both the points x i and their probabilities p i , i = 1, . . ., k) turns out to have an analytical expression for any k.For other distances, a closed-form solution is not available, but the x i and the p i can be iteratively computed by alternately solving two sets of equations till convergence, in a similar fashion to optimal quantization.It may happen that the solution is "distribution-free", meaning that the probability transforms F(x i ) and p i do not depend on the particular distribution function F selected; or, on the contrary, that the solution directly depends on F: in the latter case, we derived the sets of equations to be satisfied by the solution for a wide array of continuous random distributions.
We underline that this class of discretization techniques represents a valid alternative to the other extant procedures, among them, the consolidated moment-matching technique, whose applicability is however hindered by the often unattainable hypothesis of finiteness of the first integer moments.This class is also competitive if compared to optimal quantization and modifications thereof, since for the latter the algorithm leading to the best discrete approximation is very similar, being based on the minimization of the mean squared error between the two rvs.What we cannot expect from the suggested class of discrete approximations is the preservation of the first (and second) moment; this can be seen as a shortcoming, but only if one is used to deal with or synthesize random distributions in terms of moments or if the transformation of rvs one needs to approximate is a smooth function.In the financial or insurance fields, one is more familiar with quantiles (the most popular risk measure for market risk is the Value at Risk (VaR), which is nothing else than a quantile of a loss distribution over a fixed time horizon; see, e.g., McNeil, Frey and Embrechts 2005) and then adopting a discrete approximation which minimizes a distance between cdfs may be intuitively more appropriate and convenient.Moreover, many continuous distributions do not even possess the second or the first moment (just think of the Cauchy distribution) and then moment-matching and quantization techniques would fail in providing a discrete approximation, which is instead guaranteed by the proposed class (with rare exceptions arising if the Cramér distance is considered).The practical problem presented in the third section, concerning a complex function of several random variables, illustrates how the proposed procedures can overperform the standard moment matching approach.
Future research will examine other statistical distances between cdf, in particular asymmetrical ones, and derive the corresponding optimal k-point discrete approximation from their minimization.We will also examine possible extensions of this class of discretization techniques to the bivariate and more generally d-variate case.Although it can be naïvely judged to be straightforward, high dimensionality leads to non-negligible theoretical and computational issues: apart from the choice of the statistical distance to be minimized (an analogous quadratic distance between joint cdfs or the energy distance, see Rizzo and Székely 2016, can be considered) the most challenging matter is represented by the selection of all feasible d-variate supports for the discrete approximation (they are not necessarily a d-variate Cartesian product of univariate supports, and this complicates the evaluation of the statistical distance) and by the possible statistical dependence between the univariate components, which unavoidably calls for the use of copulas.Similar kinds of problems arise in the computation of principal points for multivariate random distributions (Flury 1990).
Another aspect that can be investigated is related to the construction of discrete approximations with "a priori" assigned support points x i .This can be useful when one wants to construct a discrete counterpart or analog to a continuous probability distribution, and not just a discrete approximation.Typically, the discrete counterpart to a continuous distribution supported on R (or R + ) is constructed as the discrete distribution supported on Z (or N) preserving the expression of the continuous cdf at the integer support points (Chakraborty 2015).Alternatively, a discrete counterpart can be constructed by considering the same integer support as above and minimizing one of the statistical distances examined in this work with respect to the probabilities p i , which now constitute a countable set.By considering any parametric random distribution with infinite support, we can thus generate several new count distributions that can be actually regarded as its discrete counterparts and can be employed for modeling count data.
and then obtain a secondorder equation in Q 1 whose unique feasible solution is Q 1 = 9+

Fig. 2
Fig. 2 Comparison between the true distribution function of the sum of exponential rvs and its three approximations (based on 7-point univariate discrete approximations); in the right panel, we zoom in on the right tail of the distribution

Table 1 k
-point discrete approximations of two continuous random distributions based on the minimization of the Cramér-von Mises distance; μ and σ 2 indicate expectation and variance of the approximation (

Table 6
KS distance between the true and the approximated distribution function of the sum of three independent exponential rvs

Table 7
Comparative closeness study for a hollow rectangular tube