Improving randomness characterization through Bayesian model selection.

Random number generation plays an essential role in technology with important applications in areas ranging from cryptography to Monte Carlo methods, and other probabilistic algorithms. All such applications require high-quality sources of random numbers, yet effective methods for assessing whether a source produce truly random sequences are still missing. Current methods either do not rely on a formal description of randomness (NIST test suite) on the one hand, or are inapplicable in principle (the characterization derived from the Algorithmic Theory of Information), on the other, for they require testing all the possible computer programs that could produce the sequence to be analysed. Here we present a rigorous method that overcomes these problems based on Bayesian model selection. We derive analytic expressions for a model's likelihood which is then used to compute its posterior distribution. Our method proves to be more rigorous than NIST's suite and Borel-Normality criterion and its implementation is straightforward. We applied our method to an experimental device based on the process of spontaneous parametric downconversion to confirm it behaves as a genuine quantum random number generator. As our approach relies on Bayesian inference our scheme transcends individual sequence analysis, leading to a characterization of the source itself.

Introduction

Random numbers have acquired an essential role in our daily lives because of our close relationship with communication devices and technology. There are also numerous scientific techniques and applications that rely fundamentally on our ability for generating such numbers and typically pseudo-random number generators (pRNGs) suffice for those purposes. A new alternative has been proposed by exploiting the inherently probabilistic nature of quantum mechanical systems. These Quantum Random Number Generators (QRNGs) are in principle superior to their classical counterparts and recent experiments have shown ref. 1 that they can reach the same quality as commercial pRNGs. However, the natural question of how to assess whether a sequence is truly random is not yet fully established. Pragmatically, the NIST test suite[2] has become the standard method for analysing sequences coming from a RNG. The suite is based on testing certain features of random sequences that are hard to reproduce algorithmically, such as its power spectrum, longest string of consecutive 1’s, and so on. Even though it constitutes an easily applicable procedure, recent findings show that its reliance on P-values is a drawback[3, 4], while its lack of formality is a major disadvantage. On the other hand, although no definition of randomness is deemed absolute, a rigorous characterization is presented by the Algorithmic Theory of Information (ATI) but it is unfortunately inapplicable in real cases[5]. An alternative which overcomes both formal and applicability issues is the Borel-normality criterion[6] (BN). Intuitively, this approach works by successively compressing a given dataset, e.g. of M bits, by taking strings of β consecutive bits and computing the frequency of occurrences of each of those possible strings. For example, β = 1 corresponds to looking for the frequencies of the strings {0, 1} in the dataset , while β = 2 corresponds to analysing the frequencies of the strings {00, 01, 10, 11}, and so on. The whole sequence is said to be Borel-normal if the frequencies are bounded individually according toand with β an integer ranging from 1 to β max = log2 log2  M. It is important to mention that BN criterion is a (nearly) necessary condition for a sequence to be considered random[5]. Note that this test is restricted to a-single-sequence classification, so it cannot determine the random character of the generating source.

In the present work, we show that randomness characterization can also be addressed using a Bayesian inference approach for model selection[7], borrowing the compression scheme of BN. For simplicity, for a fixed β we denote each string with its decimal base representation . The first step consists in identifying the models which could have generated a compressed dataset . For instance if β = 1, we can describe it as M realizations of a Bernoulli process, leading to two possible models: with and without bias. Similarly, for β = 2, a model represents a way of constructing with bias in some of the 22 possible strings. A simple combinatorial counting reveals that all the possible bias assignations correspond to all partitions of the four strings of .

Thus, in general, given the set , let denote the family of its possible partitions[8], with the Bell’s numbers and the Stirling numbers of the second kind, which counts the different ways of grouping 2 elements into K sets. Formally, would refer to the -th partition into K subsets, but for notational simplicity we will omit henceforth the index . To each partition α ( there corresponds a unique model which assigns a probability p to string according to the following rule:This means that all strings contained in a given subset ω ( are deemed equiprobable within the specified model. Thus, keeping β fixed, the likelihood of observing the given dataset in a model is:where is the frequency of string and we have defined as the aggregate frequencies of the strings in the subset ω (. (For further use, we also introduce the relative aggregate frequencies ). From this perspective, only the model that is symmetric under any reordering of the possible strings is identified with a complete random source, because any other model entails biases assignations according to the strings’ grouping represented by the corresponding partition. This symmetry only exists when the partition is the set itself, hence we denote .

Consider now that when characterising randomness the only essential feature is whether bias for or against some strings is present, but the degree of bias is irrelevant. We can eliminate the dependence on the bias parameters by multiplying with a prior for and derive the so called evidence for a given model[9]. Following[10], we use the Jeffreys prior for it yields a model’s probability distribution invariant under reparametrization and provides a measure of a model’s complexity, thus giving a mathematical representation of Occam’s Razor principle[10-12]. After integrating in the parameter space, we arrive at (see Supplementary Information (SI), Sec. 2)Eq. (4) is our main result, for it will let us perform the model selection straightforwardly. For , its evidence is fairly intuitive:Finally, we want to infer the model that best describes our source, after a dataset is given. Using Bayes’ theorem the posterior distribution reads:Henceforth we will consider a uniform prior over models (which is justified in SI), so the model’s posterior is simply proportional to its evidence.

Suppose now we want to assess whether a source can be considered truly random. This is performed in two steps. As the first step, we need a model ranking procedure based on the posterior distribution. The second step consists in quantifying the goodness of our choice of model.

As a decision rule for the ranking process we use the Bayes Factor[13] perspective,Thus, we will choose over whenever BF > 1. It has been shown that BF provides a measure of goodness of fit and if is the true model[14].

To implement the second step, which is nothing more than a hypothesis testing problem, we have two alternatives: either we check whether log10 BF ≥ 2 which is considered decisive in favour of model [13], or we compute the ratio between the posterior and the prior of a given model to assess how certain the posterior has become under the information provided by the dataset.

From a computational point of view notice that the evaluation of the posterior requires to being able to compute the normalization factor that appears in (6). When the number of models is very large we can choose either to work with a subspace of models or use the logarithm of the Bayes Factor, as in this case the normalisation factor cancels out.

It is clear that a full test of randomness requires different values of β to be used for the same dataset, while the strings should be short enough so that the M bits allow for each of the possible models to be sampled at least once. Thus, heuristically, whence we can reproduce the BN limit[6], β max ~ log2 log2 (M), after using an asymptotic expansion for the Bell number.

Note that by fixing β we have the set of parameters , whose space can be divided into regions identifying the likeliest model according to Eq. (4). As illustrative cases, in Fig. 1 we show a phase-type diagram for β = 1 and β = 2 (upper and lower panel, respectively), where the orange-filled area delimits the parameters values that renders the likeliest model. The top panel includes the bounds according to the BN criterion (green curves) given by Eq. (1), and shows that for any sequence length, M, our method allows for considerably smaller variations of γ 0. This is a significant improvement, since only necessary criteria exist for testing randomness. The lower panel depicts the analogous regions when β = 2, for which there are fifteen models (see a list in the SI) and we have fixed two frequencies: γ 1 = 1/6 and γ 2 = 1/4. The complete models distribution can be deduced from the structure of this graph, by distinguishing, a posteriori, the equiprobable strings for which the corresponding model is the likeliest. Thus more information than complete randomness classification can be readily obtained from our method.

Figure 1

Phase diagram of Randomness Characterisation. Division of the parameter space into regions according to the likeliest model. The top figure corresponds to β = 1 in terms of the frequency γ 0 of the string 0 and the sample size M. The green curves corresponds to Borel’s normality criterion, while the red curves are Borel-type bounds obtained by an approximation obtained from Eq. (4) (see Sec. 3 of SI). The bottom plot corresponds to β = 2 where each coloured area identifies the likeliest model in that region. Here we fixed the frequencies γ 1 = 1/6 and γ 2 = 1/4 and varied the frequency γ 0 of the string 00 and the sample size M.

Also in Fig. 1, the red curves of the β = 1 case are bounds obtained by comparing the likelihood of with models involving partitions into K = 2 subsets. Agreement with the regions boundary is excellent. Our choice of K = 2 is justified as we would expect that models corresponding to partitions into two subsets to be the closest ones to the model . An explicit expression for these bounds is derived in SI, Sec. 3, and Extended Data Figs 2 and 3 depict that they also bound considerably well the region in which is the likeliest for β = 2.

For further benchmarking, we have compared our method against the NIST test suite[2]. The result is depicted in Fig. 2, as a function of the sequence length M and bias b employed to generate a 0. The upper panel on Fig. 2 shows the averaged number of tests passed when employing the NIST suite, while the lower one shows the frequency of being the likeliest, for β = 1, 2 and 3. We believe that our technique can contribute to test the quality of RNG in a more stringent form, since by applying a single test thrice (once for each value of β), we determined more precisely the random character of the sample of sequences.

Figure 2

Comparison with NIST Suite test. Comparison of the bias allowed on a given sequence for it to be considered random using the NIST suite (upper panel) and our Bayesian method for randomness characterisation (lower panel).

As an application, we have tested our method in a bit sequence obtained experimentally from the differences in time detection in the process of spontaneous parametric down conversion (SPDC). Sequences generated via a SPDC photon-pair source have been shown to fulfil with ease the BN criterion, and to pass comfortably the NIST’s suite[1]. In the SPDC process a laser pump beam illuminates a crystal with a χ (2) nonlinearity, leading to the annihilation of pump photons and the emission of photon pairs, typically referred to as signal and idler[15]. Our experimental setup is shown in Extended Fig. 1 and we explain how to construct a 0 or 1 symbol from the detection signals in Section 1 of SI. We generated a 4 × 109 bits sequence, so β max ~ 4. When 1 ≤ β ≤ 3, we used all the possible models in the comparison, while, for computational ease, when β = 4, we restricted the model space to the 32, 768 models corresponding to K = 1 and K = 2 subsets (consider that ). Our inference showed that was the likeliest model for every value of β.

As explained above, to achieve a full characterization of our QRNG as a random source, we need to go further from the model ranking based on the Bayes Factor and measure our certainty that is the true model governing the source. This (un)certainty quantification is the hallmark of Bayesian statistics, since represents the probability that modelling our QRNG as a random source is correct. Computing this posterior distribution directly from Bayes’ Theorem, Eq. 6, we arrive at the values shown in Table 1 for each β. The first three values are at least 0.95, but the corresponding to β = 4 is about 0.32, considerably smaller. However, this represents an improvement of order 104 when compared with the initial value for the prior, . Alternatively, we computed log10 BFsym, for each value of β. The values reported in Table 1 correspond to the comparison of and the second likeliest model, hence the inequality for β > 2. These two criteria combined lead us to conclude that there is decisive evidence for our hypothesis that is the underlying model driving our source, thus verifying that the photonic RNG is strictly random in the sense described in the article.

Table 1

Posterior calculated for a dataset of 4 × 109 bits.

β	P(ℳ_sym|sˆ)	log10 BFsym,α′
1	0.99993	4.15
2	0.99927	≥3.55
3	0.95374	≥1.84
4	0.31862	≥3.16

From a more general perspective, we propose that quantifies our certainty on the hypothesis that a sequence was generated using the biases on strings associated with α (. Because Bayesian methods entails a model’s generalizability[9, 10], the likeliest model provides a characterization of the source of . All partitions can be identified with standard computational packages, although it can be computationally demanding for sequences of ~1010 bits. In any case, once a partition is given, its model’s likelihood is easily found using Eq. (4). A simplified analysis can be performed with the BN-type bounds given in Section 3 of the SI, which also leads to more stringent criteria than other approaches.

Supplementary information