Monte Carlo Statistical Tests for Identity of Theoretical and Empirical Distributions of Experimental Data

Often experimental work requires analysis of many datasets derived in a similar way. For each dataset it is possible to find a specific theoretical distribution that describes best the sam‐ ple. A basic assumption in this type of work is that if the mechanism (experiment) to generate the samples is the same, then the distribution type that describes the datasets will also be the same [1]. In that case, the difference between the sets will be captured not through changing the type of the distribution, but through changes in its parameters. There are some advantag‐ es in finding whether a type of theoretical distribution that fits several datasets exists. At first, it improves the fit because the assumptions concerning the mechanism underlying the experi‐ ment can be verified against several datasets. Secondly, it is possible to investigate how the variation of the input parameters influences the parameters of the theoretical distribution. In some experiments it might be proven that the differences in the input conditions lead to quali‐ tative change of the fitted distributions (i.e. change of the type of the distribution). In other cases the variation of the input conditions may lead only to quantitative changes in the output (i.e. changes in the parameters of the distribution). Then it is of importance to investigate the statistical significance of the quantitative differences, i.e. to compare the statistical difference of the distribution parameters. In some cases it may not be possible to find a single type of dis‐ tribution that fits all datasets. A possible option in these cases is to construct empirical distri‐ butions according to known techniques [2], and investigate whether the differences are statistically significant. In any case, proving that the observed difference between theoretical, or between empirical distributions, are not statistically significant allows merging datasets and operating on larger amount of data, which is a prerequisite for higher precision of the statistical results. This task is similar to testing for stability in regression analysis [3].


Introduction
Often experimental work requires analysis of many datasets derived in a similar way. For each dataset it is possible to find a specific theoretical distribution that describes best the sample. A basic assumption in this type of work is that if the mechanism (experiment) to generate the samples is the same, then the distribution type that describes the datasets will also be the same [1]. In that case, the difference between the sets will be captured not through changing the type of the distribution, but through changes in its parameters. There are some advantages in finding whether a type of theoretical distribution that fits several datasets exists. At first, it improves the fit because the assumptions concerning the mechanism underlying the experiment can be verified against several datasets. Secondly, it is possible to investigate how the variation of the input parameters influences the parameters of the theoretical distribution. In some experiments it might be proven that the differences in the input conditions lead to qualitative change of the fitted distributions (i.e. change of the type of the distribution). In other cases the variation of the input conditions may lead only to quantitative changes in the output (i.e. changes in the parameters of the distribution). Then it is of importance to investigate the statistical significance of the quantitative differences, i.e. to compare the statistical difference of the distribution parameters. In some cases it may not be possible to find a single type of distribution that fits all datasets. A possible option in these cases is to construct empirical distributions according to known techniques [2], and investigate whether the differences are statistically significant. In any case, proving that the observed difference between theoretical, or between empirical distributions, are not statistically significant allows merging datasets and operating on larger amount of data, which is a prerequisite for higher precision of the statistical results. This task is similar to testing for stability in regression analysis [3].
Formulating three separate tasks, this chapter solves the problem of identifying an appropriate distribution type that fits several one-dimensional (1-D) datasets and testing the statistical significance of the observed differences in the empirical and in the fitted distributions for each pair of samples. The first task (Task 1) aims at identifying a type of 1-D theoretical distribution that fits best the samples in several datasets by altering its parameters. The second task (Task 2) is to test the statistical significance of the difference between two empirical distributions of a pair of 1-D datasets. The third task (Task 3) is to test the statistical significance of the difference between two fitted distributions of the same type over two arbitrary datasets. Task 2 can be performed independently of the existence of a theoretical distribution fit valid for all samples. Therefore, comparing and eventually merging pairs of samples will always be possible. This task requires comparing two independent discontinuous (stair-case) empirical cumulative distribution functions (CDF). It is a standard problem and the approach here is based on a symmetric variant of the Kolmogorov-Smirnov test [4] called the Kuiper two-sample test, which essentially performs an estimate of the closeness of a pair of independent stair-case CDFs by finding the maximum positive and the maximum negative deviation between the two [5]. The distribution of the test statistics is known and the p value of the test can be readily estimated. Tasks 1 and 3 introduce the novel elements of this chapter. Task 1 searches for a type of theoretical distribution (out of an enumerated list of distributions) which fits best multiple datasets by varying its specific parameter values. The performance of a distribution fit is assessed through four criteria, namely the Akaike Information Criterion (AIC) [6], the Bayesian Information Criterion (BIC) [7], the average and the minimal p value of a distribution fit to all datasets. Since the datasets contain random measurements, the values of the parameters for each acquired fit in Task 1 are random, too. That is why it is necessary to check whether the differences are statistically significant, for each pair of datasets. If not, then both theoretical fits are identical and the samples may be merged. In Task 1 the distribution of the Kuiper statistic cannot be calculated in a closed form, because the problem is to compare an empirical distribution with its own fit and the independence is violated. A distribution of the Kuiper statistic in Task 3 cannot be estimated in close form either, because here one has to compare two analytical distributions, but not two stair-case CDFs. For that reason the distributions of the Kuiper statistic in Tasks 1 and 3 are constructed via a Monte Carlo simulation procedures, which in Tasks 1 is based on Bootstrap [8].
The described approach is illustrated with practical applications for the characterization of the fibrin structure in natural and experimental thrombi evaluated with scanning electron microscopy (SEM).

Theoretical setup
The approach considers N 1-D datasets . ≤ x n i i ) of a given random quantity under equal conditions. The datasets contain samples of the same random quantity, but under slightly different conditions.
The procedure assumes that M types of 1-D theoretical distributions are analyzed. Each of them has a probability density function PDF j (x, p → j ), a cumulative distribution function CDF j (x, p → j ), and an inverse cumulative distribution function invCDF j (P, p → j ), for j=1, 2, …,

M.
Each of these functions depends on n j p -dimensional parameter vectors p → j (for j=1, 2, …, M), dependent on the type of theoretical distribution.

Task 1 -Theoretical solution
The empirical cumulative distribution function CDF e i (.) is initially linearly approximated over (n i +1) nodes as (n i -1) internal nodes CDF e and In that rare case the upper external node will never be with positive abscissa because ( It is convenient to introduce "before-first" ples. When for some k=1,2,…,n i and for p>1 it is true that then the initial approximation of CDF e i (.) contains a vertical segment of p nodes. In that case the p nodes on that segment are replaced by a single node in the middle of the vertical segment CDF e i (x k i )=(k + p / 2 − 1 / 2) / n i . The described two-step procedure [2] results in a strictly increasing function CDF e i (.) in the closed interval The non-zero part of the empirical density PDF e i (.) is determined in the closed interval i as a histogram with bins of equal area (each bin has equal product of density and span of data). The number of bins b i is selected as the minimal from the Scott [9], Sturges [10] and Freedman-Diaconis [11] suggestions / iqr e i ). In the last three formulae, fl(y) stands for the greatest whole number less or equal to y. The lower and upper margins of the k-th bin m d ,k i and m u,k i are determined as quantiles (k-1) . The described procedure [2] results in a histogram, where the relative error of the worst PDF e i (.) estimate is minimal from all possible splitting of the samples into b i bins. This is so because the PDF estimate of a bin is found as the probability that the random variable would have a value in that bin divided to the bin's width. This probability is estimated as the relative frequency to have a data point in that bin at the given data set. The closer to zero that frequency is the worse it has been estimated. That is why the worst PDF estimate is at the bin that contains the least number of data points. Since for the proposed distribution each bin contains equal number of data points, any other division to the same number of bins would result in having a bin with less data points. Hence, the relative error of its PDF estimate would be worse. If the samples are distributed with density PDF j (x, p → j ), then the likelihood of the dataset χ i The numerical characteristics of the j th theoretical distribution fitted to the dataset χ i are calculated as: Theory and Applications of Monte Carlo Simulations • inter-quartile range: iqr j The quality of the fit can be assessed using a statistical hypothesis test. The null hypothe- , which means that the fit is not good. The Kuiper statistic V j i [12] is a suitable measure for the goodness-of-fit of the theoretical cumulative distribution functions The theoretical Kuiper's distribution is derived just for the case of two independent staircase distributions, but not for continuous distribution fitted to the data of another [5]. That is why the distribution of V from (1), if H 0 is true, should be estimated by a Monte Carlo procedure. The main idea is that if the dataset χ i = (x 1 i , x 2 i , ..., x n i i ) is distributed in compliance with the 1-D theoretical distributions of type j, then its PDF would be very close to its estimate PDF j (x, p → j i ), and so each synthetic dataset generated from PDF j (x, p → j i ) would produce Kuiper statistics according to (1), which would be close to zero [1].
The algorithm of the proposed procedure is the following: 1. Construct the empirical cumulative distribution function CDF e i (x) describing the data in χ i .

2.
Find the MLE of the parameters for the distributions of type j fitting χ i as 6. The p-value P value, j fit ,i of the statistical test (the probability to reject a true hypothesis H 0 that the j th type theoretical distribution fits well to the samples in dataset χ i ) is estimated as the frequency of generating synthetic Kuiper statistic greater than the actual Kuiper statistic V j i from step 4: In fact, (2) is the sum of the indicator function of the crisp set, defined as all synthetic datasets with a Kuiper statistic greater than V j i .
The performance of each theoretical distribution should be assessed according to its goodness-of-fit measures to the N datasets simultaneously. If a given theoretical distribution cannot be fitted even to one of the datasets, then that theoretical distribution has to be discarded from further consideration. The other theoretical distributions have to be ranked according to their ability to describe all datasets. One basic and three auxiliary criteria are useful in the required ranking.
The basic criterion is the minimal p-value of the theoretical distribution fits to the N datasets: minP value, j fit = min{P value, j fit ,1 , P value, j fit ,2 , ..., P value, j fit ,N }, for j=1, 2, ...,M . (3) The first auxiliary criterion is the average of the p-values of the theoretical distribution fits to the N datasets: The second and the third auxiliary criteria are the AIC-Akaike Information Criterion [6] and the BIC-Bayesian Information Criterion [7], which corrects the negative log-likelihoods with the number of the assessed parameters: That solves the problem for selecting the best theoretical distribution type for fitting the samples in the N datasets.

Task 2 -Theoretical solution
The second problem is the estimation of the statistical significance of the difference between two datasets. It is equivalent to calculating the p-value of a statistical hypothesis test, where the null hypothesis H 0 is that the samples of χ i1 and χ i2 are drawn from the same underlying continuous population, and the alternative hypothesis H 1 is that the samples of χ i1 and χ i2 are drawn from different underlying continuous populations. The two-sample asymptotic Kuiper test is designed exactly for that problem, because χ i1 and χ i2 are independently drawn datasets. That is why "staircase" empirical cumulative distribution functions [13] are built from the two datasets χ i1 and χ i2 : The "staircase" empirical CDF sce i (.) is a discontinuous version of the already defined empirical CDF e i (.). The Kuiper statistic V i1,i2 [12] is a measure for the closeness of the two 'staircase' empirical cumulative distribution functions CDF sce i1 (.) and CDF sce i2 (.): The distribution of the test statics V i1,i2 is known and the p-value of the two tail statistical test with null hypothesis H 0 , that the samples in χ i1 and in χ i2 result in the same 'staircase' empirical cumulative distribution functions is estimated as a series [5] according to formulae (9) and (10).
The algorithm for the theoretical solution of Task 2 is straightforward: 1. Construct the "staircase" empirical cumulative distribution function describing the data in χ i1 as CDF sce 2. Construct the "staircase" empirical cumulative distribution function describing the data 3. Calculate the actual Kuiper statistic V i1,i2 according to (8).
4. The p-value of the statistical test (the probability to reject a true null hypothesis H 0 ) is estimated as: where λ = 1 V i 1,i 2 ( n i 1 n i 2 n i 1 + n i 2 + 0.155 + 0.24 If P value,e i1,i2 <0.05 the hypothesis H 0 is rejected.

Task 3 -Theoretical solution
The last problem is to test the statistical significance of the difference between two fitted distributions of the same type. This type most often would be the best type of theoretical distribution, which was identified in the first problem, but the test is valid for any type. The problem is equivalent to calculating the p-value of statistical hypothesis test, where the null ) fitted to the datasets χ i1 and χ i2 are identical, and the alternative hypothesis The test statistic again is the Kuiper one V j i1,i2 : As it has already been mentioned the theoretical Kuiper's distribution is derived just for the case of two independent staircase distributions, but not for the case of two independent continuous cumulative distribution functions. That is why the distribution of V from (11), if H 0 is true, should be estimated by a Monte Carlo procedure. The main idea is that if H 0 is true, ), fitted to the merged dataset χ i1+i2 formed by merging the samples of χ i1 and χ i2 [1].
The algorithm of the proposed procedure is the following: 1. Find the MLE of the parameters for the distributions of type j fitting χ i1 as 3. Find the MLE of the parameters for the distributions of type j fitting χ i2 as

7.
Find the MLE of the parameters for the distributions of type j fitting χ i1+i2 as

Software
A platform of program functions, written in MATLAB environment, is created to execute the statistical procedures from the previous section. At present the platform allows users to test the fit of 11 types of distributions on the datasets. A description of the parameters and PDF of the embodied distribution types is given in Table 1 [14,15]. The platform also permits the user to add optional types of distribution.
The platform contains several main program functions. The function set_distribution contains the information about the 11 distributions, particularly their names, and the links to the functions that operate with the selected distribution type. Also, the function permits the inclusion of new distribution type. In that case, the necessary information the user must provide as input is the procedures to find the CDF, PDF, the maximum likelihood measure, the negative loglikelihood, the mean and variance and the methods of generating random arrays from the given distribution type. The function also determines the screen output for each type of distribution.

Generalized extreme value distribution Weibull distribution
Parameters μ ∈ ( − ∞; +∞), σ ∈ (0; +∞), ξ ∈ ( − ∞; +∞) Parameters x ∈ ( − ∞; +∞) (ξ = 0) Support x ∈ 0; +∞) Another key function is fitdata. It constructs the fit of each theoretical distribution over each dataset, evaluates the quality of the fits, and gives their parameters. It also checks whether two distributions of one type fitted to two different arbitrary datasets are identical. In other words, this function is associated with Task 1 and 3 (see sections 2.1 and 2.2). To execute the Kuiper test the function calls kutest. Finally, the program function plot_print_data provides the on-screen results from the statistical analysis and plots figures containing the pair of distributions that are analyzed. The developed software is available free of charge upon request from the authors provided proper citation is done in subsequent publications.

Source of experimental data for analysis
The statistical procedures and the program platform introduced in this chapter are implemented in an example focusing on the morphometric evaluation of the effects of thrombin concentration on fibrin structure. Fibrin is a biopolymer formed from the blood-borne fibrinogen by an enzyme (thrombin) activated in the damaged tissue at sites of blood vessel wall injury to prevent bleeding. Following regeneration of the integrity of the blood vessel wall, the fibrin gel is dissolved to restore normal blood flow, but the efficiency of the dissolution strongly depends on the structure of the fibrin clots. The purpose of the evaluation is to establish any differences in the density of the branching points of the fibrin network related to the activity of the clotting enzyme (thrombin), the concentration of which is expected to vary in a broad range under physiological conditions.  Fig. 1). Some main parameters of the 12 collected datasets are given in Table 2.
An automated procedure is elaborated in MATLAB environment (embodied into the program function find_distance.m) to measure lengths of fibrin strands (i.e. sections between two branching points in the fibrin network) from the SEM images. The procedure takes the file name of the fibrin image (see Fig. 1) and the planned number of measurements as input.
Each file contains the fibrin image with legend at the bottom part, which gives the scale, the time the image was taken, etc.
The first step requires setting of the scale. A prompt appears, asking the user to type the numerical value of the length of the scale in μm. Then the image appears on screen and a red line has to be moved and resized to fit the scale (Fig. 2a and 2b). The third step requires a red rectangle to be placed over the actual image of the fibrin for selection of the region of interest (Fig. 2c). With this, the preparations of the image are done, and the user can start taking the desired number of measurements for the distances between adjacent nodes (Fig. 2d).
Using this approach 12 datasets containing measurements of lengths between branching points of fibrin have been collected (Table 2) and the three statistical tasks described above are executed over these datasets.

Task 1 -Finding a common distribution fit
A total of 11 types of distributions (  (3) and (4), and the highest of (5) and (6), i.e. it is the worst fit. The types of distributions worth using are mostly the lognormal distribution (having the lowest AIC and BIC), and the generalized extreme value (having the highest possible meanP value, j fit ). Figure 3 presents 4 of the 11 distribution fits to DS4. Similar graphical output is generated for all other datasets and for all distribution types. , and the minimal values for AIC j and BIC j across the datasets for each distribution type. The bold and the italic values are the best one and the worst one achieved for a given criterion, respectively.   Table 4 contains the p-value calculated according to (9) for all pairs of distributions. The bolded values indicate the pairs, where the null hypothesis fails to be rejected and it is possible to assume that those datasets are drawn from the same general population. The results show that it is possible to merge the following datasets: 1) DS1, DS2, DS3, D4 and DS8; 2) DS7, DS10, and DS11; 3) DS5 and DS12. All other combinations (except DS5 and DS10) are not allowed and may give misleading results in a further statistical analysis, since the samples are not drawn from the same general population. Figure 4a presents  =0.2830). The Kuiper statistic for identity of the empirical distributions, calculated according to (8), is V 4,9 =0.5005, whereas according to (9) P value,e 4,9 =2.024e-24<0.05. Therefore the null hypothesis is rejected, which is also evident from the graphical output. In the same fashion, Figure 4b presents (8), is V 1,4 =0.1242, whereas according to (9) P value,e 1,4 =0.1957>0.05. Therefore the null hypothesis fails to be rejected, which is also confirmed by the graphical output.

Task 3 -Identity of fitted distributions
As concluded in task 1, the lognormal distribution provides possibly the best fit to the 12 datasets. Table 5 contains the p-values calculated according to (12) for the lognormal distribution fitted to the datasets with 1000 Monte Carlo simulation cycles. The bold values indicate the pairs, where the null hypothesis fails to be rejected and it is possible to assume that the distribution fits are identical. The results show that the lognormal fits to the following datasets are identical: 1) DS1, DS2, DS3, and DS4; 2) DS1, DS4, and DS8; 3) DS7, DS10, and DS11; 4) DS5 and DS10; 5) DS5 and DS12. These results correlate with the identity of the empirical distribution. Figure 5a presents the fitted lognormal distribution over =0.5120, mode 7 9 =0.4322, std 7 9 =0.2442, iqr 7 9 =0.2926). The Kuiper statistic for identity of the fits, calculated according to (11), is V 7 4,9 =0.4671, whereas according to (12), P value, 7 4,9 =0<0.05. Therefore the null hypothesis is rejected, which is also evident from the graphical output. In the same fashion, Fig. 5b (11), is V 7 1,4 =0.03288, whereas according to (12), P value,7 1,4 =0.5580>0.05. Therefore the null hypothesis fails to be rejected, which is also evident from the graphical output.  The statistical procedures described above have been successfully applied for the solution of important medical problems [16; 17]. At first we could prove the role of mechanical forces in the organization of the final architecture of the fibrin network. Our ex vivo exploration of the ultrastructure of fibrin at different locations of surgically removed thrombi evidenced gross differences in the fiber diameter and pore area of the fibrin network resulting from shear forces acting in circulation (Fig. 6). In vitro fibrin structures were also generated and their equivalence with the in vivo fibrin architecture was proven using the distribution analysis described in this chapter (Fig. 7). Stretching changed the arrangement of the fibers (Fig. 7A) to a pattern similar to the one observed on the surface of thrombi (Fig. 6A); both the median fiber diameter and the pore area of the fibrins decreased 2-3-fold and the distribution of these morphometric parameters became more homogeneous (Fig. 7B). Thus, following this verification of the experimental model ultrastructure, the in vitro fibrin clots could be used for the convenient evaluation of these structures with respect to their chemical stability and resistance to enzymatic degradation [16]. B. Fiber diameter (upper graphs) and fibrin pore area (lower graphs) were measured from the SEM images of the DG thrombus shown in A using the algorithms described in this chapter. The graphs present the probability density function (PDF) of the empirical distribution (black histogram) and the fitted theoretical distribution (grey curves). The numbers under the location of the observed fibrin structure show the median, as well as the bottom and the top quartile values (in brackets) of the fitted theoretical distributions (lognormal for fiber diameter and generalized extreme value for area). The figure is reproduced from Ref. [16].

Figure 7.
Changes in fibrin network structure caused by mechanical stretching. A. SEM images of fibrin clots fixed with glutaraldehyde before stretching or following 2-and 3-fold stretching as indicated, scale bar = 2 μm. B. Fiber diameter (upper graphs) and fibrin pore area (lower graphs) were measured from the SEM images illustrated in A using the algorithms described in this chapter. The graphs present the probability density function (PDF) of the empiric distribution (black histogram) and the fitted theoretical distribution (grey curves). The numbers under the fibrin type show the median, as well as the bottom and the top quartile values (in brackets) of the fitted theoretical distributions (lognormal for fiber diameter and generalized extreme value for area). The figure is reproduced from Ref. [16].
Application of the described distribution analysis allowed identification of the effect of red blood cells (RBCs) on the structure of fibrin [17]. The presence of RBCs at the time of fibrin formation causes a decrease in the fiber diameter (Fig. 8) based on a specific interaction between fibrinogen and a cell surface receptor. The specificity of this effect could be proven partially by the sensitivity of the changes in the distribution parameters to the presence of a drug (eptifibatide) that blocks the RBC receptor for fibrinogen (compare the median and interquartile range values for the experimental fibrins in the presence and absence of the drug illustrated in Fig. 8). It is noteworthy that the type of distribution was not changed by the drug, only its parameters were modified. This example underscores the applicability of the designed procedure for testing of statistical hypotheses in situations when subtle quantitative biological and pharmacological effects are at issue.

Discussion and conclusions
This chapter addressed the problem of identifying a single type of theoretical distribution that fits to different datasets by altering its parameters. The identification of such type of distribution is a prerequisite for comparing the results, performing interpolation and extrapolation over the data, and studying the dependence between the input parameters (e.g. initial conditions of an experiment) and the distribution parameters. Additionally, the procedures included hypothesis tests addressing the identity of empirical (stair-case) and of fitted distributions. In the case of empirical distributions, the failure to reject the null hypothesis proves that samples come from one and the same general population. In the case of fitted distributions, the failure to reject the null hypothesis proves that although parameters are random (as the fits are also based on random data), the differences are not statistically significant. The implementation of the procedures is facilitated by the creation of a platform in MAT-LAB that executes the necessary calculation and evaluation procedures.
Some parts of the three problems analyzed in this chapter may be solved using similar methods or software tools different from the MATLAB procedures described in section 3. Some software packages solve the task of choosing the best distribution type to fit the data [18,19]. The appropriateness of the fit is defined by the goodness-of-fit metrics, which may be selected by the user. The Kolmogorov-Smirnov statistics is recommended for the case of samples with continuous variables, but strictly speaking the analytical Kolmogorov-Smirnov distribution should not be used to calculate the p-value in case any of the parameters has been calculated on the basis of the sample as explicitly stated in [19]. Its widespread application, however, is based on the fact that it is the most conservative, i.e. the probability to reject the null hypothesis is lower compared to the other goodness-of-fit criteria. Some available tools [20] also use analytical expressions to calculate the p-value of the Kolmogorov-Smirnov test in the case of a sample that is normally distributed, exponentially distributed or extreme-value distributed [21,22]. Those formulae are applied in the lillietest MATLAB function from the Statistical toolbox, where Monte-Carlo simulation is conducted for the other distributions. It is recommended to use Monte-Carlo simulation even for the three aforementioned distributions in case any of the parameters has been derived from the sample. Some applications calculate a goodness-of-fit metrics of a single sample as a Kuiper statistics (e.g. in the awkwardly spelled kupiertest MATLAB function of [23]) and the p-value is calculated analytically. The main drawback of that program is that the user must guarantee that the parameters of the theoretical distribution have not been calculated from the sample. Other available applications offer single-sample Kuiper test (e.g. v.test function in [24]) or single-and two-sample Kuiper tests (e.g. KuiperTest function in [25]), which use Monte-Carlo simulation. The results of the functions v.test and KuiperTest are quite similar to those presented in this chapter, the main difference being our better approximation of the empirical distribution with a linear function, rather than with a histogram. Our approach to calculate p-values with Monte-Carlo simulation stems from the previously recognized fact that "…if one or more parameters have to be estimated, the standard tables for the Kuiper test are no longer valid …" [26]. Similar concepts have been proposed by others too [27].
An advantage of the method applied by us is that the Kuiper statistics is very sensitive to discrepancies at the tails of the distribution, unlike the Kolmogorov-Smirnov statistics, whereas at the same time it does not need to distribute the data into bins, as it is for the chi-square statistics. Another advantage is that the method is very suitable for circular probability distributions [23,24], because it is invariant to the starting point where cyclic variations are observed in the sample. In addition it is easily generalized for multi-dimensional cases [25].
A limitation of our method is that it cannot be used for discrete variables [25], whereas the Kolmogorov-Smirnov test could be easily modified for the discrete case. The second drawback is that if the data are not i.i.d. (independent and identically distributed), then all Bootstrap and Monte-Carlo simulations give wrong results. In that case, the null hypothesis is rejected even if true, but this is an issue with all Monte-Carlo approaches. Some graphical and analytical possibilities to test the i.i.d. assumption are described in [19].
Further extension of the statistical procedures proposed in this chapter may focus on the inclusion of additional statistical tests evaluating the quality of the fits and the identity of the distributions. The simulation procedures in Task 3 may be modified to use Bootstrap, because this method relies on fewer assumptions about the underlying process and the associated measurement error [28]. Other theoretical distribution types could also be included in the program platform, especially those that can interpret different behaviour of the data around the mean and at the tails. Finally, further research could focus on new areas (e.g. economics, finance, management, other natural sciences) to implement the described procedures.