Deterministic sampling from uniform distributions with Sierpiński space-filling curves

In this paper the problem of sampling from uniform probability distributions is approached by means of space-filling curves, a topological concept that has found a number of important applications in recent years. Departing from the theoretical fact that they are surjective but not necessarily injective, the investigation focused upon the structure of the distributions obtained when their domains are swept in a uniform and discrete manner, and the corresponding values used to build histograms, that are approximations of their true PDFs. This work concentrates on the real interval [0,1] and the Sierpiński space-filling curve was chosen because of its favorable computational properties. In order to validate the results, the Kullback–Leibler and other divergence measures are used when comparing the obtained distributions in several levels of granularity with other already established sampling methods. In truth, the generation of uniform random numbers is a deterministic simulation of randomness using numerical operations. In this fashion, sequences resulting from this sort of process are not truly random. Despite this, and to be coherent with the literature, the expression “random number” will be used along the text to mean “pseudo-random number”.


Introduction
In the last decades the proliferation of low cost and faster digital computers has expanded the development and applicability of techniques using stochastic simulation.The uncertainty in such systems is usually simulated by means of random numbers.Although some time ago the purpose of the simulations was aimed at qualitative understanding, the explosive availability of processing power in the last decades made it possible to reach more quantitative results.In this fashion, large scale simulations made it feasible the estimation of objects with large variances and very precise estimation of variables under study.Accordingly, random numbers with high quality and long periods are needed-portability and efficiency are essential as well.Frequently, simulations of probabilistic models need random variables with specific probability distributions and, in practice, many algorithms for generation of these non-uniform random variables are based on certain transformations of uniform random numbers.Therefore, uniform random numbers can be considered the basis for probabilistic simulation.In truth, uniform random number generation is a deterministic process which imitates randomness by means of arithmetic operations.Therefore, sequences obtained in this way are radically different from sequences of truly random numbers.The notion of absolute randomness is an idealized concept, and there were initiatives in the sense of "statistically purify" it by means of models of physical processes.While there is some logical basis for that, this kind of approach is not very well accepted in the field of simulation.This is so because numbers generated by using a particular physical process are not truly random numbers, which exist only conceptually.As another example we have the argument that there is no ideal dice without any bias, capable of simulating Bernoulli trials.In addition, there is the infeasibility of reproducing the same sequence of random numbers several times and in different contexts.This type of procedure is fundamental for checking simulation results obtained by others, and code debugging.
However, generating true randomness on digital computers is not a feasible task, taking into account their deterministic nature.In this fashion, sequences obtained from deterministic processes are known as pseudorandom ones.There are at least two different ways of dealing with this situation: the pragmatic one, that is, to use pseudorandom numbers with some optimism or just to give up usinig randomness on computer simulations and converging to the concept of derandomization.The first alternative will demand the execution of statistical tests on every set produced by the chosen generators in different applications.So, they should be "stressed" in several directions, so that randomness be at least approximated, even though never reached.In Marsaglia (1993) it is stated that all deterministic methods for producing randomness fail in some application, and only experience and imagination of developers and users can lead to a better understanding of this type of event.Despite all those true statements, considerable progress has been made since the 1940s, when computer-generated random numbers were successfully used in the implementation of the so-called Monte Carlo methods.
For the sake of efficiency, typical algorithms have been based on recurrence relations, which can be viewed as discrete dynamical systems "equipped" with finite memory.These two characteristics of discreteness and finiteness show up whenever considering any type of method for digital computer generation of random objectsactually, present day computers are unable to exactly represent irrational numbers, leaving available only the rational ones, up to a limited numerical precision.Finiteness, by definition, implies in the production of periodic sequences.For example, sometimes it is possible to find people with the belief that certain nonlinear discrete dynamical systems are candidates for deterministic schemes of random number generation.But, on digital computers, it is simply not possible to simulate chaotic behavior at all, considering the finite precision of numbers processed by present day machines, among other things.So, the main problem in random number generation on computers is quick generation, using small amounts of memory, and producing sequences with huge periods.
An alternative approach to handle the difficult situation described above, by means of pseudorandom numbers, is to look for more rational and precise mathematical solutions.According to von Neumann (1963), "It is true that a problem that we suspect of being solvable by random methods may be solvable by some rigorously defined sequence."In the same direction, some researchers proposed alternative methods using deterministic versions of Monte Carlo methods, called quasi-Monte Carlo methods.Also, and in the same line of reasoning, this article presents a deterministic and practical method for asymptotically sampling from a 1-dimensional uniform probability distribution with arbitrary precision levels.The proposal is based on the space-filling curve synthesized by Sierpiński (Sagan 1994), coupled to some simple transformations described below.
Many applications of space-filling curves and their approximations have been recently proposed in the literature (Goertzel 1999;Lera and Sergeyev 2010;Sergeyev et al. 2013).Probably their most important and amazing property is the ability to completely fill compact regions of higher dimensional spaces, in the sense of (typically) being surjective and continuous at the same time.

Sierpi ński space-filling curves
It is possible to define space-filling curves as surjective and continuous functions from [0,1] to compact subsets of finite dimensional vector spaces, usually identified to R n .These objects were well-studied in the past, with many theoretical results establishing necessary conditions for their existence (Sagan 1994).In addition, several important mathematicians defined concrete examples and established several interesting properties (Sagan 1994) of SFCs.Decades later, researchers have found significant applications of space-filling curves to several fields, including global optimization of numerical functions and design of certain devices (Goertzel 1999).However, implementation using digital computers causes certain difficulties,mainly due to their finite word length.Space-filling curves are often defined by means of infinite expansions and use, for instance, the property that elements in [0,1] can be represented in the form 0.t 1 t 2 t 3 t 4 ... in a given basis B, where each t i is an integer between 0 and B-1.Accordingly it is necessary to find precise approximations of SFCs to pursue this kind of approach in real world computers.The most adequate candidate seems to be the Sierpiński curve, taking into account the availability of its precise and simple defining formulas, as shown below (Sagan 1994; Oliveira and Petraglia 2011): where f is a bounded, even and continuous real function whose expression is given by The resulting curve (x(t), y(t)) is a 2-dimensional SFC showing good behavior in numerical computations.When constructing higher-dimensional SFCs, certain results are fundamental.Before stating them, however, it is necessary to present some definitions.
Definition 1 A function ϕ : [0, 1] → R is uniformly distributed with respect to the Lebesgue measure if, for any (Lebesgue) measurable set A ⊂ R, we have where Λ 1 is the Lebesgue measure in the real line.
Considering that the Sierpiński SFC is measure-preserving ( Sagan 1994, page 111), it is clear that its coordinate functions are uniformly distributed and stochastically independent, and can be used to obtain higher-dimensional SFCs with coordinates that are not constant, being continuous and stochastically independent.It occurs that, for higher-dimensional sets, approximations of original Sierpiński based SFCs were not able to completely fill up the associated compact domains.This is due to distortions caused by numerical approximations, despite the theoretical curve being a space-filling one.
In order to find a more robust curve, it is possible to compose the Sierpiński function with an invertible measure-preserving transformation that is a natural extension of a particular generalized Lüroth series transformation (Dajani and Kraaikamp 2002), mapping [0,1)x[0,1) onto itself.Finally, a new SFC mapping and vice-versa.By denoting such a transformation by T , its definition is It is worth highlighting that despite Steinhaus' theorem is stated only for continuous functions, it is also valid for surjective (over [0,1]), piecewise continuous coordinate functions.In this fashion, the following result is true.

Proposed method
Considering that space-filling curves are, by definition, surjective, it is obvious that every single point in their codomains is "visited" at least once.The fundamental question is: what kind of distribution corresponds to their full excursion in the respective domain?
In principle, each type of SFC has one specific solution for this issue, and the corresponding theoretical treatment is not easy, for sure.Turning back to the problem at hand, and focusing on the unidimensional case, let us consider a Sierpiński SFC having  Marsaglia and Tsang (2000) In Tables 1, 2, 3 and Figs. 1, 2 and 3, numerical results corresponding to all divergences are displayed, evidencing the convergence of PDFs to the uniform distribution -1000 bins were used for the construction of histograms.Also, according to simulations it is possible to observe that the SFC-based distribuition is faster in the convergence process.
In addition, for the sake of showing more detailed information, histogram sets displayed in Figs. 4, 5, 6, 7 through 8 are shown in the sequence, resulting from sampling according to the chosen methods.The sequence of figures clearly shows the convergence process taking place, and closes the pictorial presentation.Notice that histograms corresponding to SFCs are "smoother" than the others.In these simulations, only 50 bins were used, for the sake of better visualization.Finally, the graph displayed in Fig. 9 corresponds to the PACF (Partial Autocorrelation Function) of a sequence formed by 1000000 samples, with 60000 lags.As can be seen, the picture presents a typical configuration of sequences of pseudorandom numbers, with no detectable cyclic or recursive behavior.

Conclusions
Through several simulations, this paper has shown that it is possible to (approximately and) deterministically sample from U (0, 1) by means of Sierpiński space-filling curves.The proposed method was compared to 4 other well-established algorithms and the "proximity" to the uniform distribution was measured with the Kullback-Leibler,  Hellinger and Bhattacharyya divergences, which are widely accepted as good indices of discrepancy between PDFs.By analyzing the numerical results, it is possible to infer that the presented algorithm converges faster than the other methods.This type of result may be very useful in applied fields, like cryptography, for example.

Fig. 2
Fig. 2 Hellinger divergences to uniform distribution for each method

Fig. 8
Fig. 3 Bhattacharyya divergences to uniform distribution for each method