Interval-based, nonparametric approach for resampling of fuzzy numbers

In this paper, we propose two new nonparametric resampling methods for the simulation of bootstrap-like samples of fuzzy numbers. The generated secondary samples are based on an input set (i.e., a primary sample) consisting of left–right fuzzy numbers. The proposed approaches utilize random simulations in a way which, to some extent, resembles a bootstrap. However, contrary to the classical bootstrap approach, the proposed methods are based on alpha-cuts of fuzzy numbers, which are generated in a new nonparametric way. Therefore, these procedures give us an opportunity to create ”not exactly the same as previous” fuzzy numbers and also lead to greater diversity of the obtained output. Moreover, we check whether the introduced methods can be successfully applied in two statistical tests about the mean value of a population of fuzzy numbers.


Introduction
Nowadays, computer-aided simulations constitute an important tool for solving practical problems in many areas, like, physics, mathematics, biology, chemistry. However, in order to simulate a model of a real process that has a random component, we need algorithms, which generate random variables in a specified way. Two important approaches should, in particular, be mentioned: the Monte Carlo methods and resampling algorithms. In the first case, a sample of iid (independent, identically distributed) random variables is used to solve the given problem, which is too complex for analytical evaluation (see, e.g., Robert and Casella (2004) for additional details and further discussion). In the second case, an initial sample of observed values is reused to create a secondary sample in order to, e.g., calculate the standard error of Communicated  The John Paul II Catholic University of Lublin, ul. Konstantynów 1H, 20-708 Lublin, Poland some complicated test statistics (see, e.g., Efron (1982) for an introduction and various examples). Resampling techniques, commonly known as bootstrap methods, are especially useful in inferential problems of mathematical statistics. They are used for estimation of probability distributions of sample statistics when analytical methods are either too complex or unavailable.
The same applies to the problems in which we have to deal with complex phenomena that are characterized not only by randomness, but also by fuzzy imprecision, as well. In such cases, for the description of the considered models we may use the concept of a fuzzy random variable that may be regarded as an imprecise (fuzzy) counterpart of a well known crisp (e.g., real) random variable.
Simulation methods for fuzzy random variables strongly depend upon their interpretation. There exist several mathematical models of fuzzy random variables. For their description, the reader is advised to read pertaining literature, such as, e.g., a very good monograph by Couso et al. (2014) or overview papers (Gil and Hryniewicz 2009;Gil et al. 2006a). The most popular interpretation of fuzzy random variables, known as "epistemic", is based on the model proposed in the papers by Kwakernaak (1978Kwakernaak ( , 1979. In this model, a fuzzy random variable describes imprecise (fuzzy) perception of an unobserved crisp random variable. The "epistemic" model of fuzzy random variables has been applied for solving many real-life problems using computer-aided simulations. Its applications were described in numerous papers. For example, the authors of this paper used it to solve problems from such different areas as: pricing of financial and insurance instruments Romaniuk 2013, 2017), estimation of the maintenance costs of a water distribution system (Romaniuk 2016(Romaniuk , 2018 or Bayesian statistical decisions in reliability . The simulations which were considered in these papers, usually consisted in the generation of random hidden crisp origins, and respective membership functions (e.g., in the form of triangles with edges of random length). Successful applications of simulation methods in the case of "epistemic" fuzzy random variables can be explained using the results by Hryniewicz (2015), who noticed that fuzzy random variables, defined according to the definition by Kwakernaak, can be described in a fully probabilistic way using infinitely dimensional probability distributions.
The second popular definition of a fuzzy random variable was proposed in the seminal paper by Puri and Ralescu (1986). This definition is based on the notion of set-valued mapping and random sets, and its interpretation is called "ontic." Simulation of "ontic" fuzzy random variables is much more difficult than in the case of "epistemic" ones. The main reason for this is the nonexistence of the concept of a probability distribution that describes fuzzy random "ontic" observations. We may define classical probability distribution only for certain sample characteristics, such as, the sample mean, but not for fuzzy observations. Moreover, for "ontic" random fuzzy variables popular measures of variability (such as variance) do not exist, and other measures should be used instead. All these interpretational problems make simulation processes of fuzzy random variables having "ontic" interpretation much more difficult. For example, Colubi et al. (2002) considered simulation methods for different types of both one-and multidimensional fuzzy variables in this setting. They used these methods for the analysis of asymptotic behavior of a fuzzy arithmetic mean, expressed in terms of the strong law of large numbers, and of the law of iterated logarithm. The process of simulation itself was thoroughly examined in the paper by González-Rodríguez et al. (2009). They proposed two different approaches, based on the concept of support functions. The first one is related to simulations of Hilbert space-valued random elements with a projection on the cone of all fuzzy sets. The second one imitates the representation of elements of a separable Hilbert space for an orthonormal basis directly on the space of fuzzy sets. Both of these approaches were compared, and their comparison showed that the second method is more adequate for modeling realistic situations.
The lack of a natural probability distribution of "ontic" fuzzy random variables makes resampling (bootstrap) methods a valuable alternative to Monte Carlo simulation. Indeed, bootstrap methods have been successfully used in statistical tests about the expected value of a fuzzy random variable (see, e.g., Gil et al. 2006b;González-Rodríguez et al. 2006;Montenegro et al. 2004), and in other types of statistical tests in fuzzy environment (see, e.g., Ramos-Guajardo et al. 2010;Ramos-Guajardo and Lubiano 2012). In the aforementioned papers, bootstrap samples enable the authors of the considered statistical tests to estimate a nominal significance level of a test via an empirical percentage of rejections of a true null hypothesis. In this approach, a bootstrap-based estimator serves as an empirical benchmark for the considered statistical test. Another bootstrap method, namely the weighted bootstrap, was used by Hung (2006) in the construction of the minimum inaccuracy fuzzy estimator, the calculation of its standard error, and the construction of appropriate confidence intervals.
A classical bootstrap approach has one disadvantage, appearing when the original fuzzy sample is small. In such a case, a bootstrap sample consists of few distinct values. This could be considered as an obstacle when this sample has to be used in modeling of complex phenomena, described, e.g., by complex functions of fuzzy random variables. Therefore, a modification of this method, with the aim of increasing the diversity of simulated results, seems to be needed. This is the main goal of this paper, in which we propose a new approach for the simulation of quasi-bootstrap populations of fuzzy random variables. Our approach is applicable to both "epistemic" and "ontic" fuzzy random variables.
Our idea to simulate more diverse fuzzy random populations is implemented in this paper in the form of two new nonparametric resampling methods for the so-called left-right fuzzy numbers (LFRNs). The proposed approach consists of two steps. In the first step, we use classical resampling methods in which a primary sample of fuzzy observations is reused by us in order to randomly generate the secondary sample. In the second step, we make some random perturbations of membership functions of fuzzy elements of the secondary sample. Therefore, our algorithms generate fuzzy numbers that may differ from the fuzzy numbers included in the original primary sample. Our new methods of resampling may be considered as a kind of bootstrap-like generation methods. The fuzzy numbers generated in the secondary sample may be thus "not exactly the same", but rather "similar" (in some mathematical sense) to the values from the primary sample. Therefore, they have greater diversity, when compared to the secondary samples simply resampled from primary samples, but still possess the same main statistical characteristics.
In order to measure the aforementioned similarity and the diversity of a new secondary sample, we compare the generated sample and the input set using four types of measures of similarity and two types of triangular fuzzy numbers. The applicability of the strong law of large numbers and the law of iterated algorithm as indicators of convergence of the generated samples is also checked. Of course, LRFNs generated using proposed methods should have some practical value and should be applicable in solving the real-life problems. Therefore, we discuss an application of the introduced methods in two bootstrapped versions of statistical tests about the mean of a population of fuzzy numbers. An empirical pvalue of these tests serves as a benchmark in the performed comparisons.
The paper is organized as follows. In Sect. 2, basic definitions of fuzzy sets and random fuzzy variables have been recalled. Moreover, we have presented the descriptions of statistical tests used for testing the hypotheses about the expected value. Next, in Sect. 3, we describe the proposed algorithms for the generation of bootstrap-like secondary samples. Then, in Sect. 4 we describe the results of the experimental verification of the properties of the proposed procedures. The application of the proposed new bootstrap procedures in statistical testing has been presented in Sect. 5. The paper is concluded in its last section.

Fuzzy numbers and random fuzzy numbers
Let us present basic definitions and notation, concerning the simulation of fuzzy random variables, which will be used in this paper. Additional details can be found in, e.g., Gil and Hryniewicz (2009) and Gil et al. (2006a).
Definition 1 A fuzzy numberã is a fuzzy subset of R for which μã is a normal, upper-semicontinuous, fuzzy convex function with a compact support.
Then,Ã(0) is the closure of the set {x : μÃ (x) > 0}. A fuzzy numberã is a fuzzy subset of R for which μã is a normal, upper-semicontinuous, fuzzy convex function with a compact support. Then, for each α ∈ [0, 1], the α-level set a(α) is a closed interval of the formã A left-right fuzzy number (which is further abbreviated as LRFN) is a fuzzy number with the membership function of the form A triangular fuzzy numberã, denoted further by a L , a C , a R , is an LRFN with the membership function of the form where a L is the left end of its support, a C -its core, and a R -the right end of its support. Some examples of triangular fuzzy numbers are shown in Fig. 2. Fuzzy random variables are generalizations of ordinary (crisp) random variables. Historically, the first widely accepted definition of the fuzzy random variable was proposed by Kwakernaak (1978Kwakernaak ( , 1979. Below, we present this definition in the version elaborated by Kruse and Meyer (1987). Kruse and Meyer 1987) Let ( , A, P) be a probability space, where is the set of all possible outcomes of the random experiment, A is a σ -field of subsets of (the set of all possible events of interest), and P is a probability measure associated with ( , A). A mapping X :

Definition 2 (See
is the space of all fuzzy numbers, is called a fuzzy random variable if it satisfies the following properties: , are usual real-valued random variables associated with ( , A, P).
The fuzzy random variable defined according to Definition 2 has a seemingly natural interpretation. It may be considered as a fuzzy perception of an unknown true realvalued random variable associated with a random experiment and referred to as 'the original' of the considered fuzzy random variable. This interpretation is called "epistemic" and allows to process fuzzy data, interpreted as realizations of fuzzy "epistemic" random variables, in relatively easy way. For example, "epistemic" interpretation of fuzzy random variable allows us to simulate fuzzy random data in a relatively direct way, i.e., without making many additional assumptions.
Another, and from the mathematical point of view more general, definition was proposed by Puri and Ralescu (1986). Puri and Ralescu 1986). Given a probability space ( , A, P), a mapping X :

Definition 3 (See
→ F c (R) is said to be a fuzzy random variable (also referred to as random fuzzy set) if for each α ∈ [0, 1] the set-valued mapping X α : → K c (R), where K c (R) is the class of the non-empty compact intervals and X α (ω) = (X (ω)) α for all ω ∈ , is a compact convex random set (that is, a Borel-measurable mapping with respect to the Borel σ -field generated by the topology associated with the Hausdorff metric on K c (R)).
The fuzzy random variable defined by Definition 3 may be used for the analysis of random events (random data) that are intrinsically fuzzy. This interpretation is called "ontic" and may be used to process random data presented in the form of fuzzy random sets. Such fuzzy data exist in practice (see, e.g., the book by Viertl (2011) for examples), but their analysis is much more difficult. Because of these difficulties, there exist problems with the simulation of fuzzy random "ontic" data.

Measures of similarity
To compare some properties of two fuzzy numbers, like shape or location of their membership functions, one can use various measures of similarity. In this paper, we apply three classical measures: the supremum, the l 1 metric, the Hausdorff metric for fuzzy sets (see, e.g., Zwick et al. (1987) for additional details), and a more complex distance measure, which was introduced by Tran and Duckstein (2002). Then, all of these measures will be used to compare the LRFNs generated using the methods proposed in Sect. 3 with the fuzzy numbers taken from an initial (primary) sample.
Ifã andb are fuzzy sets, then the supremum similarity measure, introduced in Nowakowska (1977), is defined for their membership functions μã(x) and μb(x), as In the case of the l 1 metric, proposed in Kaufman (1975), an appropriate measure is given by There are various ways to extend the Hausdorff distance to a metric for fuzzy sets (see, e.g., Zwick et al. 1987). In this paper, we will use this distance in the version proposed by Ralescu and Ralescu (1984), and defined as The fourth distance measure considered in this paper was introduced by Tran and Duckstein (2002), and it is given by the following formula where w(α) is a certain weighting function. In this paper, we assume that w(α) = 1, so each α-cut in the measure (1) has the same significance (see Tran and Duckstein (2002) for other possible types of the weighting function and further discussion). Of course, many other types of measures of similarity between two fuzzy sets have been proposed, and some of them are widely used. But, in the following, we will focus our attention only on the previously mentioned four measures of similarity. Three of them (namely m ∞ , m l 1 , m H ) represent very "classical" and "standard" approaches. They have strict relationships with mathematical measures of similarity, which are known and used for crisp values (i.e., real numbers), or with geometrical measures applied to points in space. Therefore, they are precisely and intuitively understood and easy to implement. Characteristics of many other measures are directly compared to the properties of these three ones. Obviously, these measures have some disadvantages, too. For example, the measure m ∞ takes into account only the supremum of the difference between the membership functions of two fuzzy sets. Therefore, even two "very similar" (in a broad sense of this word) fuzzy sets can significantly differ, if this measure is applied.
We also consider the fourth measure, denoted by m TD . According to its proposers (Tran and Duckstein 2002), this measure is specifically tailored for the LRFNs, which are our main objects of interest in this paper. Moreover, in Tran and Duckstein (2002), the authors enumerated many advantages of this measure, such as straightforward computation, facility of interpretation for decision makers, robustness, flexibility, and transitivity (see also Zhu and Lee 1992). After some comparisons, Tran and Duckstein conclude in Tran and Duckstein (2002) that this measure is at least as reasonable as other existing ones. From our point of view, the application of m TD in this paper seems to be a "one step further" (i.e., beyond "very classical" approaches) in comparison with the usage of the three previously mentioned measures.
In Sect. 4, we apply these four measures to the check of similarity and diversity of triangular fuzzy numbers, generated using the methods introduced in Sect. 3 and using the classical bootstrap approach. Of course, other types of measures can be also utilized for this purpose. However, quite surprisingly, the obtained results seem to behave in a very stable way without unexpected differences, regardless of the measure used. Therefore, it seems to us that application of other measures of similarity should lead to more or less the same overall results.

Tests of the fuzzy mean value
There are many types of statistical tests for an expected value of a fuzzy random variable (see, e.g., Gil et al. 2006b;González-Rodríguez et al. 2006;Körner 2000;Montenegro et al. 2004). We focus on only two of them, which will be used in Sect. 5 as examples of application of the introduced nonparametric simulation methods.
The first considered test is an asymptotic test introduced in Körner (2000). Let us assume thatã is an LRFN with a core, which is given by a single value. Then, we have m a = a L (1) = a R (1), l a = m a − a L (0), r a = a R (0) − m a .
The d 2 distance between two LRFNsã andb is defined as The values of R 1 and R 2 are defined analogously (see Körner 2000). For this type of distance, we have the following corollary, which was proved in Körner (2000): Corollary 1 Let X 1 , X 2 , . . . , X n be a sample of LRFNs. Then where ξ 1 , ξ 2 , ξ 3 are independent N (0, 1)-distributed random variables and λ 1 , λ 2 , λ 3 are the eigenvalues of the matrix Moreover, an asymptotic test of the hypothesis where w 2 q is the q-th quantile of an ω 2 distribution with respect to the eigenvalues λ 1 , λ 2 , λ 3 .
The above-mentioned ω 2 distribution has a rather complex structure, which is known only for some special cases (see Körner 2000).
The second considered test was developed in González-Rodríguez et al. (2006) and Montenegro et al. (2004). It is based on a metric introduced in Bertoluzza et al. (1995), which was generalized in Körner and Näther (2002). The D ϕ W metric for two LRFNsã,b is defined as where , and W , ϕ are two weighting normalized measures (see Bertoluzza et al. (1995) for some examples of W , ϕ and further details). Then, we have the following corollary, which was established in González-Rodríguez et al. (2006): . . , X n be a sample of LRFNs. In testing, the null hypothesis where z q is the q-th empirical quantile of the bootstrap distribution, which is given by and with where X * 1 , X * 2 , . . . , X * n is a bootstrap sample obtained from the initial sample X 1 , X 2 , . . . , X n .
We will use both these tests in the experimental analysis of the practical applicability of our simulation procedures.

Generation of the secondary (bootstrap) sample
Let A = {ã 1 , . . . ,ã m } be a primary sample of LRFNs. These values are treated as an input set for the methods proposed further on in this paper. We assume that we do not have (and, moreover, we do not need) any additional information about a source (or a model) of the fuzzy numbers belonging to A. Note, however, that in many cases, which are known from literature, such additional information is often assumed (see, e.g., Colubi et al. 2002;Hryniewicz 2015;Hryniewicz et al. 2015;Romaniuk 2013, 2017;Romaniuk 2016 for various approaches to the problem of fuzzy numbers modeling). Therefore, only a strictly nonparametric way should be used to build a secondary (bootstrap) sample B = {b 1 , . . . ,b n } of LRFNs, which should be, in some way, "similar" to the fuzzy numbers from A. Letã j (α) = a L j (α), a R j (α) be an α-cut ofã j for some α ∈ [0, 1]. For simplicity, we assume that there are k + 1 Algorithm 1: Initialization procedure Input: A primary sample A = {ã 1 , . . . ,ã m }, the number of possible α-cuts k + 1. Output: Sets of the cores and spreads for A. Find a set of values of the cores C(1) = {a 1 (1), . . . , a m (1)} and order it increasingly; However, this requirement can be easily relaxed in a simulation procedure presented further.
During the first step of an initialization procedure (a setup of simulation, see Algorithm 1), a set of cores C(1) is found, based on A. Hence, we have For simplicity of notation, we assume that the set C(1) is already ordered, i.e., a 1 (1) ≤ a 2 (1) ≤ · · · ≤ a k (1).
During the second step of the initialization procedure, sets of incremental spreads for all possible α-cuts are constructed. Let be the difference between left ends of α-cuts for α i+1 and α i , for the given fuzzy numberã j . We call such a difference an incremental left spread for the level i. In the same manner, we have which is the difference between the right ends of α-cuts for α i and α i+1 , for the given fuzzy numberã j . It will be called an incremental right spread for the level i. Then, the sets of left and right incremental spreads, given by for α k−1 , α k−2 , . . . , α 0 can be found. It should be noted that the construction of (5) has to be made from the highest value of α to the lowest one (i.e., from the core of a fuzzy number to its support). We also assume, in the same manner as for the set of cores C(1), that each of the sets (5) is already ordered, so that for all α i . Let us illustrate this initialization procedure with a numerical toy-example.

The d-method based on a discrete distribution d(x)
Let us start from the description of a generation procedure in case of the d-method, based on a discrete probability distribution d(x). In the proposed procedure, two steps are necessary to construct a fuzzy numberb j ∈ B, where j = 1, . . . , n (see also Algorithm 2). Firstly, the value of a core b j (1) is found, using a uniform discrete distribution for the values from the set C(1). It means that the generated value b j (1) = C is a random element, taken from the set C(1), according to the probability distribution d(x). In this paper, we assume that d(x) is Algorithm 2: Secondary sample generation for the d(x)method Input: Sets of the cores and the incremental spreads for A, the number of LRFNs in a secondary sample n, the number of possible α-cuts k + 1. Output: A secondary sample B generated using the discrete distribution d(x). for j ← 1 to n do Randomly draw a value of a core b j (1) from the set C(1), using the discrete uniform distribution for m elements; Constructb j from the obtained α-cuts and append it to B; end uniform on C(1), i.e., where l = 1, . . . , m. Therefore, we randomly (and uniformly) pick up a single value from the set C(1) and treat it as the core of the new, constructed LRFNb j . Secondly, consecutive α-cuts of the givenb j are found, starting from its core and ending at its support. Thus, we proceed fromb j (α k−1 ) down tob j (0). For each α i , the value of the left end of the α-cut ofb j is found, using where S L (α i ) is an independently drawn random value from the set S L (α i ). Once again, the uniform discrete distribution d(x) is used, for which where l = 1, . . . , m. In the same manner, the right end of each α-cut ofb j is constructed, using where S R (α i ) is independently drawn from the set S R (α i ), using the same uniform discrete distribution d(x). Formulas (6) and (7) mean that the new left (or right, respectively) end of α i -cut is constructed, based on subtracting (or adding) a random element from the set S L (α i ) (or S R (α i )) from (to) the previously generated left (right) end of the α i+1 -cut. Therefore, this new fuzzy numberb j is approximated using intervals for the consecutive values of α (from 1 at the top to 0 at the bottom).
In some sense, the fuzzy numberb j , obtained in this way, is similar to the LRFNs from the primary sample A. The core ofb j is one of the "true" cores from C(1), and its spreads are drawn from the "true" spreads belonging to S L (α i ) or S R (α i ). It is easily seen that we have (1), so the expected values of the core and the spreads ofb j are precisely equal to the respective means for LRFNs from A.
In the same way, meaning thatb j exactly "imitates" the statistical behavior of the samples from A, without the necessity of introducing any additional knowledge about the model, which (perhaps) generates the primary sample. Now, let us continue our example by showing how the secondary bootstrap-like sample is constructed. We will show the construction of only one element of this sample. The remaining elements are constructed in the same way.
Example 1 (Continued) The core of a new element of the secondary sample B is, in this example, randomly chosen (with equal probabilities 1/3) from the set {1, 2.5, 3.5}, and let this chosen value be equal to b L 1 (1) = b R 1 (1) = 1. Then, we take randomly (also with equal probabilities 1/3) the left and right incremental spreads on the remaining α-level. Suppose, that for α = 0 we have chosen S L 1 (0) = 1.5, S R 1 (0) = 2.5. Thus, the respective α-cuts of the new elementb 1 of the secondary sample, calculated according to (6)-(7), are defined by the following limits: Fig. 4). It appears thatã 1 is "the most similar" tob 1 . Then, using the measures considered in

The w-method based on a mixed discrete uniform distribution w(x)
We can be also interested in an additional level of "freedom" when creating the secondary sample B. It is easy to see that ifb j is generated using the method described in Sect. 3.1 (i.e., using the uniform discrete distribution d(x)), its core is exactly equal to one of the values from C(1). Also its spreads are given by the respective values from the sets Yet, in some cases, creation of a more diversified sample B can be useful. Due to such diversification, the values from B could be "closer" to the (unknown) hidden model, than the samples from A, especially if the number of elements in A is strictly limited. Consider, for example, the case when there are only two fuzzy numbers in A, described by only two α-cuts. The random numbersb j , generated using the method described in Sect. 3.1, have no more than two possible values of a core and four possible left / right ends of its support. Moreover, if a more classical resampling method is taken into account (like the "classical" bootstrap), then these two elements from A are repeated infinitely often during the construction of the LRFNs from B. Thus, no new "knowledge" about other possible outcomes, which could be possibly "produced" by the unknown model, can be obtained.
Of course, notwithstanding the introduction of the diversification, the secondary sample B should be still enough "similar" to the primary set A. If such a requirement is not fulfilled, then our knowledge resulting from B can be misleading, and our suppositions about the original source (i.e., the model of A) can be incorrect. But no strict prior knowledge about the model for the primary sample was previously assumed in this paper. Therefore, the proposed generation method should be strictly nonparametric, without any additional and more detailed assumptions.
In statistics, if we do not want to introduce any prior knowledge, we have to use the so-called non-informative probability distributions. A commonly used model of such a distribution is a uniform density for an interval [c, d], denoted further by U ([c, d]). We will use this density in the construction of the probability distribution used for generation purposes.

The w(x) distribution, and its properties
Let be a strictly increasing sequence of m values, without their repetitions. We propose a density w(x), which is the composition of a discrete random distribution, and a continuous probability density, given by the following formula where and δ x (.) is the Dirac measure. If X ∼ w(x), where w(x) is given by (9), then X = x 1 or X = x m is taken with an atomic probability 1 2m . Hence, the first value x 1 or the last one x m from the sequence (8) is selected with equal probabilities. Otherwise, one of the intervals [x l−1 , x l ], for l = 1, . . . , m − 1, is designated with atomic probability 1 m . When such a single interval is selected, we have X ∼ w l−1,l (x), so the output x is generated using the uniform density U ([x l−1 , x l ]), which is described by (10).
Therefore, w(x) can be seen as a certain generalization of the discrete distribution, discussed in Sect. 3.1. The pdf w(x) also generates values from the same interval [x 1 , x m ], but they are more diversified-apart from the values directly equal to the ones from the sequence (8), all x ∈ [x 1 , x m ] can now be obtained.
Statistical characterizations of the density w(x) are summarized in the following lemma: Lemma 1 Let X ∼ w(x), where w(x) is a pdf described by (9) and (10). Then and Var X = 1 m Proof From (9) and (10), we have which concludes the proof.
From Lemma 1 we see that if X ∼ w(x), then the expected value of X is precisely equal to its meanx. But the variance s 2 w of X is not equal to the classical estimator, i.e., the standard sample variance s 2 . The difference between the variances s 2 w and s 2 can be important for the intended diversity of the LRFNs in the second sample B. We have which leads to the following remark: Remark 1 If m → ∞ and x i > 0, then s 2 w − s 2 ≥ 0. Therefore, the variability (measured by variance) of X ∼ w(x) is not lesser than the variability of X ∼ d(x), if only the size of a sample is large enough, and all x i > 0.

Generation procedure
Now, to generate a fuzzy numberb j ∈ B, if j = 1, . . . , n, instead of the discrete distribution d(x), the previously introduced density w(x) is used (see also Algorithm 3). However, an overall procedure of the construction ofb j is similar to the previous case, which is described in Sect. 3.1.
During the first step, the value of a core b j (1) = C is drawn, using the distribution w(x), based on the elements from the set C(1), so C ∼ w(x), where x ∈ C(1). Next, consecutive α-cuts of the givenb j are calculated, starting Algorithm 3: Secondary sample generation using the w(x)-method Input: Sets of the cores and the incremental spreads for A, the number of LRFNs in a secondary sample n, the number of possible α-cuts k + 1. Output: A secondary sample B generated using the distribution w(x). for j ← 1 to n do Randomly draw a value of a core b j (1) from the set C(1), using the density (9) for m elements; for i ← k − 1 to 0 do Randomly draw a value of a left incremental spread S L (α i ) from the set S L (α i ), using the density (9) for m elements; Find the left end of the Randomly draw a value of a right incremental spread S R (α i ) from the set S R (α i ), using the density (9) for m elements; Find the right end of the Constructb j from the obtained α-cuts and append it to B; end from the value α k−1 and ending at α 0 = 0. For each α i , a value of the left end ofb j (α i ) is equal to (6), where S L (α i ) is an independently drawn random value from the set S L (α i ), using the distribution w(x) for the set S L (α i ). In the same way, the right end ofb j (α i ) is given by (7), where S R (α i ) is independently drawn from the set S R (α i ), using the respective distribution w(x) for this set.
Let us continue our example using the method of generation, described in this subsection. We will show the construction of only one element of this sample. The remaining elements are constructed in the same way.
Example 1 (Continued) Let us start from the generation of the core of a new fuzzy numberb 1 . According to the density function w(x) defined by (9), there are four possibilities for choosing this value: take 1 with probability 1/6, take a randomly chosen (using the uniform distribution) number from the interval [1, 2.5] with probability 1/3, take a randomly chosen number from the interval [2.5, 3.5] with probability 1/3, take 3.5 with probability 1/6 (see also Fig. 5). Suppose that the second option has been chosen, and a new core has been set to b 1 (1) = 1.75. Now, consider the left and right spreads for α = 0. For choosing the value of the left incremental spread there are also 4 possibilities: take 1 with probability 1/6, take a randomly chosen number from the interval [1, 1.5] with probability 1/3, take a randomly chosen number from the interval [1.5, 2.5] with probability 1/3, take 2.5 with probability 1/6. Suppose that the first option has been chosen, and a new left incremental spread has been set to S L 1 (0) = 1. Similarly, for choosing the value of the right incremental spread there are 4 possibili- ties: take 1.5 with probability 1/6, take a randomly chosen number from the interval [1.5, 2] with probability 1/3, take a randomly chosen number from the interval [2, 2.5] with probability 1/3, take 2.5 with probability 1/6. Suppose that the third option has been chosen, and a new right incremental spread has been set to S R 1 (0) = 2.1. Finally, the new generated element of the secondary sample is the fuzzy number defined by its core b 1 (1) = 1.75, and its α-cut is defined by the respective limits b L 1 (0) = 0.75, b R 1 (0) = 3.85 (see Fig. 6). This new elementb 1 seems to be similar to bothã 1 andã 2 . But, using the considered measures, we can compare the pairsã 1 ,b 1 andã 2 ,b 1 , and try to decide, which of the fuzzy numbers from the initial sample is more similar tob 1 . We get m ∞ ã 1 ,b 1 = 0.75, m l 1 ã 1 ,b 1 = 1.3625, m H ã 1 ,b 1 = 0.8, m TD ã 1 ,b 1 = 1.11778 and m ∞ ã 2 ,b 1 = 0.5, m l 1 ã 2 ,b 1 = 1.29375, m H ã 2 ,b 1 = 0.95, m TD ã 2 ,b 1 = 1.23722. Then, if we choose the supremum or the TD measure, we conclude that the fuzzy numbersã 1 ,b 1 are the most similar. If, instead, we use m l 1 or m H , the triangular fuzzy numbersã 2 ,b 1 are the most similar. Let us start from a certain population P n 0 , which consists of n 0 LRFNs. From this population, we randomly draw m elements. Let these elements constitute a primary sample A m . Afterwards, using the fuzzy numbers from this primary sample, three methods (i.e., the classical bootstrap, the d-method, and the w-method) are used to generate the secondary sample B n , which consists of n elements.
In our numerical experiments, different settings are used: a moderate population P 100 (for which n 0 = 100) together with a small primary sample A 5 (where m = 5) and a moderate secondary sample B 100 (where n = 100), and a bigger population P 200 with a moderate primary sample A 100 and a rather big secondary sample B 200 . This allows us to compare outcomes for the classical bootstrap, the d-method, and the w-method, for the cases when preliminary information about a model (which is available only via the analysis of the primary sample) is very sparse (in the case of A 5 ) or relatively abundant (for A 100 ).
For simplicity, only triangular fuzzy numbers will be considered, i.e., only two α-cuts (where α 0 = 0 and α 1 = 1) are used to construct the whole LRFN. Actually, such simple types of fuzzy numbers are frequently used by practitioners. However, both the d-method and the w-method can be easily used to generate the second sample, even if more α-cuts are considered.
In the following numerical experiments, two types of triangular numbers are considered as a model for the population P n 0 . The first one (which is further referred to as the "type A number") is a fuzzy number with an expected symmetrical spread, where the center is random and has the standard normal distribution N (0, 1), and the semiwidths of the support are given as independent Chi-square variables with one degree of freedom. A similar LRFN is discussed in a detailed way in Colubi et al. (2002). The second kind (the "type B number") of a fuzzy number has a strictly non-symmetrical shape. In this case, the center points are described by the gamma distribution with the shape parameter 1 and the scale parameter 2, and the semiwidths of the support being drawn from independent exponential distributions with parameter 1 (for the left spread) or 2 (for the right spread).
We are interested in the analysis of mutual relations between the primary and secondary samples for the different generation procedures and the mentioned types of LRFNs. Therefore, properties of both the primary sample and the generated, secondary set, are statistically summarized using the sample mean versus the population mean, which are calculated for the support and the center of fuzzy numbers. From the statistical point of view, also variability of the simulated fuzzy numbers is very important. Thus, the standard deviation is also found for the support and the center of LRFN in the case of the second (i.e., generated) sample.
Moreover, the simulated fuzzy numbers should give us some additional "insight" into the model, which is (generally) completely unknown and "hidden" in the data from the primary sample. In an ideal situation, LRFNs from the second sample should be (in some way) "similar" to the numbers from the primary sample, but, simultaneously, not exactly "the same" as the elements from A m , and also "very close" to the population. Therefore, values of some measures (see Sect. 2) are evaluated for each possible pair of fuzzy numbers. These pairs consist of one "old" LRFN (i.e., from A), and one "new", generated fuzzy number (i.e., an element from B). The obtained measure values are also summarized using common tools, like minimum, maximum, mean, and standard deviation. Afterwards, we can conjecture whether some generation method produces fuzzy numbers which are "the same as", "similar" or only "close" (and to what extent) to the LRFNs from the primary sample.

Small primary sample, type A fuzzy number
Based on the small sample A 5 of the type A of triangular fuzzy numbers, three moderate secondary samples B 100 were generated, using the classical bootstrap, the d-method, and the w-method. Then, the means for the coreX * C (see Fig. 7), the left end of the supportX * L (see Fig. 8), and the right end X * R (see Fig. 9) for each of the simulated samples were calculated. From now on, the results obtained with the bootstrap are marked by circles in the graphs, with the d-method-by diamonds, and with the w-method-by squares. Horizontal bold lines correspond to the means of the primary sample A for the coreX A C , the left end of the supportX A L and the right endX A R , and axes of respective graphs start exactly in the means for the population (for the coreX C , for the left end of the supportX L and for the right endX R ).
As it is seen, each of the simulation methods behaves generally well. In each case, after generation of 30-40 fuzzy numbers, the mean of the secondary sampleX * approaches the respective mean of the primary setX A . Moreover, appli- cations of the d-method or the w-method seem to have some advantages, when they are compared to the classical bootstrap. For example, the means for these two approaches are, in general, closer toX A (i.e., the mean of A 5 ), than in the case of the bootstrap. The respective graphs are also much smoother. Surprisingly, though, in the case of the w-method, the respective mean is also closer to "the real" result-the mean of our unknown model, i.e., the population P 100 . Apart from the comparison of the means, variability of the generated LRFNs should be also considered. Hence, standard deviations for the core (see Fig. 10), the left end of the support (see Fig. 11) and its right end (see Fig. 12) are plotted. These graphs are marked in the same way as the previous ones. In each of these cases, the standard deviation of the secondary sample is the lowest when the w-method is used. Now, we compare the three secondary samples, which are generated using the considered simulation procedures, but with the help of the measures, which were recalled in Sect. 2. Let us assume that l(ã i ,b j ) is a value of some measure of similarity l(., .) between LRFNsã i ∈ A andb j ∈ B. Then, the following notation is used The respective measures of similarity are summarized in Table 1 (when B 100 is simulated using the bootstrap approach), Table 2 (in the case of the d-method) and Table 3 (for the w-method). Of course, the bootstrap only repeats fuzzy numbers, which are already present in the primary sample. Therefore, MinMin and MinMax values for the measures m l 1 , m ∞ and m H are strictly equal to zero. But in the case of d-method, even the values of these measures are more diversified, so we have StDevMin > 0. The same applies for the w-method. Therefore, these two methods produce LRFNs, which are more diversified (and "not exactly the same" in some way, too) than the numbers from A 5 . However, the generated LRFNs are also "similar" (in the sense of the applied measures of similarity) to the fuzzy numbers from the primary sample, because the obtained MinMin and MeanMin values are very close to zero. It seems that using the w-method is more promising than the d-method, because MinMax, MaxMax, and MeanMax values are generally lesser for this first approach, and MeanMin values are very similar. Hence, even LRFNs, which are "maximally" distant from the fuzzy numbers from the primary sample, are "closer" in the case of the w-method than for the d-method. Let us analyze the observed similarity in another way. In order to do this, an additional independent sample T 200 , which consists of 200 fuzzy numbers of type A, was generated. Then, three secondary sets B 200 are sampled based on A 5 , using the bootstrap, the d-method, and the w-method. We find an LRFN from each of B 200 , which is the nearest to some fuzzy number from T 200 in the sense of one of the measures m l 1 , m ∞ , m TD , m H , i.e., a value wheret i ∈ T 200 andb j ∈ B 200 , is calculated. The obtained minimal values of these measures for respective pairs of LRFNs are given in Table 4, and for each of the measure its minimum appears in boldface. As it is seen, if the wmethod is used, then the generated fuzzy number is the most similar to some element from T 200 . In some way, this new independent sample T 200 gives an additional insight into the "true model", because it is a supplementary sample from the unknown source, which models our LRFNs. Therefore, the w-method produces fuzzy numbers, which are the nearest to this model in the considered case. Note, that because the bootstrap only repeats elements from the primary sample, then for this method the obtained values of the measures are even 6-7 times bigger than for the best match.

Small primary sample, type B fuzzy number
Now we analyze the three considered simulation procedures for the case when the small primary sample A 5 consists of the strictly non-symmetrical triangular fuzzy numbers (i.e., the previously mentioned LRFNs of "type B"). The graphs of the means (for the core-see Fig. 13, for the left end of the support-see Fig. 14, for the right end of the support-see Fig. 15) are very similar to the case, which was described in Sect. 4.1. Once again, these means for the d-method and the w-method are, in general, closer to the respective means of the primary sample, than in the case of the bootstrap. Their graphs are also very smooth. Moreover, the plots of the standard deviations behave reasonably well (for the core-see Fig. 16, for the left end of the support-see Fig. 17, for the right end of the support-see Fig. 18). The obtained values are the lowest when the w-method is used. Also the characteristics of the similarity measures, which were introduced in Sect. 4.1, can be found in this case (see Tables 5, 6, 7 for the respective summaries for the different simulation approaches). In general, the conclusions are similar as in the case of type A fuzzy numbers, i.e., the bootstrap only repeats LRFNs from the primary sample, and the dmethod and the w-method produce a more diversified output, which is still similar (in the sense of the considered measures) to the values from A 5 . But the decision as to whether the d-method or the w-method is better suited, when the "maximum" distance criterion is taken into account, is not so straightforward now. As it is seen, MinMax values are lower for the d-method, but MaxMax and MeanMax are lower in the case of the w-method.  obtained minimal values of measures can be found in Table 8. Also in this case, the w-method generates fuzzy numbers, which are the most similar to some element from the set T 200 , apart from the measure m ∞ , for which the d-method gives the best result. The classical bootstrap gives values, which are even 2-3 times bigger than the best matches.

Moderate primary sample
In practical situations, apart from small statistical samples, which consist of only few values, larger samples are also used. Therefore, we also analyze the behavior of a moderate primary sample, for which m = 100 (i.e., A 100 ), and the corresponding simulated secondary sample B 200 , which is rather a big one, especially when compared to the previous examples (now we have n = 200). As it turns out, general conclusions for both type A and type B of LRFNs are very similar to the outcomes for the small sample, which were summarized in Sects. 4.1 and 4.2. Hence, we omit a more detailed discussion, in order to present another, but, in some way, supplementary approach.
Up till now, we have analyzed the speed of convergence of the mean of the secondary sampleX * to the "true" (but, in general, unknown) mean of the populationX . And, in our reasoning, three "focal points" (a core, a left and a right end of a support) have been taken into account. In Colubi et al. (2002) the authors consider an application of LIL (the law of iterated logarithm) as a tool for a convergence diagnosis for the simulated fuzzy numbers. Therefore, we will also analyze the behavior of distance between √ n/ √ 2n log log nX * andX as a function of the secondary sample size n. To keep consistency with our previous analysis, the three mentioned "focal points" will be still in the center of our attention. Hence, the distance for the core and similarly defined measures for the left and right ends of a support, will be further used, instead of the supremum distance, which is considered in Colubi et al. (2002). Because the secondary sample B 200 is rather big, the convergence speed for (11) and the rest of similar measures, as functions of n, is now more visible. We restrict our analysis only to type A fuzzy numbers, but the obtained conclusions are also similar for type B. The calculated distances as functions of the secondary sample size are plotted in Figs. 19 (the core), 20 (the left end of the support) and 21 (the right end of the support). As it is seen, the bootstrap approach is the worst one, especially for larger values of n, because the obtained distances are, in general, bigger for this simulation method. Both the d-method and the w-method produce the relatively well behaving output.

New bootstrap-like sample as a tool in statistical tests
Apart from the statistical properties of the simulated LRFNs, the possibility of applying the proposed methods in practical statistical cases was also investigated. Two types of tests for the expected value of the fuzzy numbers were considered (see Sect. 2.3 for additional details and notation) as a respective example.
The first one is a bootstrapped version of the test proposed in Körner (2000) (see Corollary 1). From now on, it will be called the K test (from its author's name) for the expected value. The second test is the procedure developed in González-Rodríguez et al. (2006) and Montenegro et al. (2004) (see Corollary 2). It will be called the GRMCG-test for the expected value (also based on the authors' names). In this case, we apply the standard uniform density as the weight normalized measure ϕ in D ϕ W (ã,b) metric (2) (see Bertoluzza et al. (1995) and Montenegro et al. (2004) for additional details and other approaches).
As an initial sample in each of these tests, three types of triangular fuzzy numbers are simulated. Types A and B are described here in Sect. 4. Type C, which was considered in Körner (2000), is a fuzzy number, in which the random center has the standard normal distribution N (0, 1), and the spreads of the support are independently drawn from the standard uniform distribution U ([0, 1]).
For each of these types of fuzzy numbers, three different simulation procedures (the classical bootstrap, the d-method and the w-method) are used to generate an input random sample for the test. The number of elements n in such a sample is varied, so that both small and medium sample sizes are considered, i.e., we set n = 5, 10, 30, 100. Also, a few values of the number of bootstrap replications r (namely r = 100, 200, 1000) are used to generate the respective bootstrapped distribution of the test statistics. In this way, we investigate the possible influence of this parameter. In each of these experiments, the whole resampling procedure is iterated 100,000 times (see, e.g., Gil et al. (2006b), González-Rodríguez et al. (2006), Montenegro et al. (2004) and Ramos-Guajardo and Lubiano (2012) for additional details of a similar approach).
Based on the respective statistics, in each of the tests of the expected value, an empirical percentage of rejectionsp at the nominal significance level p = 0.05 for the true null hypothesis is then computed. This estimated value is widely used as a benchmarking tool for the bootstrapped version of the statistical tests (see, e.g., Gil et al. 2006b;González-Rodríguez et al. 2006;Montenegro et al. 2004;Ramos-Guajardo et al. 2010;Ramos-Guajardo and Lubiano 2012). The three considered simulation procedures can then be directly compared.
In general, the simulated values ofp for all of the approaches are very close to one another, and the overall properties are very similar. In particular, the empirical percentages of rejections converge to one another for larger values of n and r (like n = 100 and r = 1000). However, there are also some significant differences. In order to emphasize them, for each experiment the value ofp, which is the nearest to the true value of the significance level p, is given in boldface.
Let us start from the K test of the expected value. As it is seen for the fuzzy numbers of type A (see Table 9), type B (see Table 10) and type C (see Table 11), a comparison of the simulation approaches seems to be quite simple. In each of these cases, the d-method leads top, which is the nearest to the assumed significance level p, apart from a few exceptions. For all of these exceptions, the classical bootstrap approach gives the most "true" answer. But even in these cases, the differences between the empirical percentages of rejections for the d-method and the classical bootstrap are not very significant (about 0.001-0.002). And these differences favor the d-method especially for smaller values of n and r . Altogether, the classical bootstrap occupies the second place with respect to the measure of proximity betweenp and p.
For the GRMCG-test, the analysis of differences between p and p is not so straightforward. In the case of type A of fuzzy numbers (see Table 12),p seems to be the nearest to the true significance level for the w-method, when a smaller number, like n = 5, 10, of the elements in the initial sample is taken into account, or for the d-method-for larger values n = 30, 100. Especially for the small samples, the classical bootstrap approach gives the worst answers and the differences between the bootstrap and one of the other approaches are quite important (about 0.008-0.01).
However, when type B fuzzy numbers are analyzed (see Table 13), the picture is not quite clear. Firstly, for n = 5, 10, the estimated percentages of rejections favor the classical bootstrap approach, because the other approaches give bigger  In the case of type C fuzzy numbers (see Table 14), it seems that the d-method or the w-method produce the most accurate estimators ofp. This can be seen especially for the smaller samples (n = 5, 10), when the classical bootstrap approach gives an estimator of the rejection rate, which is by about 0.004 smaller than for the other methods. For the largest sample (n = 100), the d-method is favored, but once again, the differences among the simulated values ofp are quite small.
Taking into account the whole analysis, it is not possible to point out the undoubtedly best simulation procedure, which gives the most accurate values ofp. However, application of the d-method or the w-method looks promising, especially for smaller initial samples.

Conclusions
In this paper, we propose two simulation algorithms for the generation of sets, consisting of LRFNs, namely the dmethod and the w-method. Both of these algorithms are based on the resampling paradigm and utilize the primary sample of fuzzy numbers in order to randomly generate the secondary bootstrap-like sample. This generation is based on α-cuts of LRFNs and a strictly nonparametric approach, without necessity of making additional assumptions about the source (or the model) of the primary sample. Our contribution in this article is fourfold. Firstly, two new numerical algorithms for the simulation of samples of the LRFNs are considered. These algorithms, similarly as the classical bootstrap methods, utilize a primary (initial) sample of random fuzzy numbers in order to generate secondary (bootstrap) fuzzy random samples. But, contrary to the classical bootstrap, these simulated secondary sets consist of values, which are "not exactly the same" as in the initial sample. In the first method, the modified direct method (called the d-method and described by a discrete probability distribution d(x)), information about the α-cuts of the LRFNs from the primary set is used. In the second method (called the w-method), a mixed discrete uniform probability distribution w(x) is used for generation purposes. In this approach, the information about the α-cuts of the observations from the primary sample is modified in a certain way, using a non-informative uniform distribution. Both of the proposed methods are used to generate the sets of LRFNs, whose diversity is in a certain sense greater than the diversity of observations from the primary sample. However, this greater diversity has been achieved without incorporation of any additional and specific assumptions about the general probability model for the initial population. Hence, both of these approaches are strictly nonparametric ones.
Secondly, the outputs for these two methods are analyzed, using the most important statistical measures. For both small and moderate primary samples, and two types of triangular fuzzy numbers, we check if the generated secondary (bootstrap-like) samples imitate well the statistical behavior of the initial population. In order to do this, the mean and the standard deviation are calculated, and applicability of the strong law of large numbers and of the law of iterated logarithm have been confirmed. We also compare the simulated secondary samples for the two proposed methods with the output of the classical bootstrap approach. It seems that the application of d(x) and w(x) distributions in bootstrapping is very promising, because the generated triangular numbers "mimic" the values from the initial sample very well. Moreover, if the previously mentioned statistical measures are taken into account, these generated values sometimes behave even better than in the case of the classical bootstrap approach applied to the same primary samples.
Thirdly, for the same sizes of samples, and two types of triangular fuzzy numbers, we check whether the simulated values are "close enough" to the fuzzy numbers from the initial set. The level of this proximity is measured using four types of measures (the supremum measure, the l 1 metric, the Hausdorff distance extended to the metric, and the measure proposed by Tran and Duckstein in Tran and Duckstein (2002)). Once again, the obtained results have been compared with the outcomes for the classical bootstrap approach. The analysis performed confirms the proposition that fuzzy numbers, generated using d(x) and w(x) distributions are very close to observations from the primary sample. Therefore, the two simulation procedures, introduced in this paper, can be used to form the secondary (bootstrap-like) sample, which is "similar", but also, in some way, different, in comparison with the initial set of observations.
And finally, we check whether these two new simulation algorithms can be successfully applied for solving some practical statistical problems. As an example, we have applied the d-method and the w-method in two statistical tests about the mean value of a population of fuzzy numbers. In these two tests, the outputs for both small and moderate primary samples have been analyzed for three types of triangular fuzzy numbers. As previously, we have compared three simulation procedures (the classical bootstrap and two methods introduced by us). In all considered cases, the difference between the nominal significance level of the test and the empirical percentage of rejections of the true null hypothesis is used as a benchmark. Once again, the algorithms introduced in this paper show their promising potential, because the difference mentioned is usually lower for the proposed bootstrap-like procedures, based on the d(x) or w(x) distributions, than for the classical bootstrap of fuzzy random variables.
The proposed methods have, in comparison with the classical bootstrap, one disadvantage, appearing when considered fuzzy numbers have their "natural" limits (e.g., when their supports must contain only nonnegative numbers). In such a case, it may happen that some of the generated elements of the secondary (bootstrap-like) sample may not fulfill such requirements. In such a case, one can introduce certain modifications of the proposed method (e.g., a simple curtailment) in order to eliminate such "unnatural" observations. However, the consequences of such modifications are yet not known and require consideration in the future research.
It should be noted that fuzzy sets introduced by Zadeh are still the most popular tool used for modeling non-random uncertainty (imprecision). There exist many extensions of fuzzy sets which can also be used for this purpose. For example, interval-valued fuzzy sets (IVFS), introduced independently by four different authors in 1975-see Nowak and Hryniewicz (2018) for references, can be used in situations, when the membership function of a fuzzy set cannot be precisely defined. Another very popular extension, widely known under the name of intuitionistic fuzzy sets (IFS), was introduced in Atanassov (1986) and can be used when we describe imprecision in terms of membership and nonmembership functions. Many of these different methods are interrelated or even formally equivalent (see, e.g., Deschrijver and Kerre 2003). Probabilistic models for IVFS and IFS variables have been already proposed in the literature. However, statistical methods for the analysis of such imprecise data practically do not exist. Notable exceptions (Akbari and Arefi 2013;Hesamian and Akbari 2017) are devoted to the analysis of IFS random data. Complex description of IVFS and IFS data makes inferential statistical procedures for these types of data extremely difficult. Therefore, simulation methods, considered in this paper, after some necessary modifications, could be applied in the statistical analysis of such data.

Compliance with ethical standards
Conflict of interest All authors of this paper declare that they have no conflict of interest.
Human participants or animals This article does not contain any studies with human participants or animals performed by any of the authors.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.