Fabrizio Biondi

Yusuke Kawamoto

Axel Legay

Louis-Marie Traonouez

Under consideration for publication in Formal Aspects of Computing Hybrid Statistical Estimation of Mutual Information and its Application to Information Flow

Keywords: Mutual information, Statistical estimation, Quantitative information flow, Hybrid method, Confidence interval, Statistical model checking

Analysis of a probabilistic system often requires to learn the joint probability distribution of its random variables. The computation of the exact distribution is usually an exhaustive precise analysis on all executions of the system. To avoid the high computational cost of such an exhaustive search, statistical analysis has been studied to efficiently obtain approximate estimates by analyzing only a small but representative subset of the system's behavior. In this paper we propose a hybrid statistical estimation method that combines precise and statistical analyses to estimate mutual information, Shannon entropy, and conditional entropy, together with their confidence intervals. We show how to combine the analyses on different components of a discrete system with different accuracy to obtain an estimate for the whole system. The new method performs weighted statistical analysis with different sample sizes over different components and dynamically finds their optimal sample sizes. Moreover, it can reduce sample sizes by using prior knowledge about systems and a new abstraction-then-sampling technique based on qualitative analysis. To apply the method to the source code of a system, we show how to decompose the code into components and to determine the analysis method for each component by overviewing the implementation of those techniques in the HyLeak tool. We demonstrate with case studies that the new method outperforms the state of the art in quantifying information leakage.

Introduction

In modeling and analyzing software and hardware systems, the statistical approach is often useful to evaluate quantitative aspects of the behaviors of the systems. In particular, probabilistic systems with complicated internal structures can be approximately and efficiently modeled and analyzed. For instance, statistical model checking has widely been used to verify quantitative properties of many kinds of probabilistic systems [START_REF] Legay | Statistical model checking: An overview[END_REF].

The statistical analysis of a probabilistic system is usually considered as a black-box testing approach in which the Fig. 1. Decomposition of a probabilistic program with input x ∈ {x 1 , x 2 , x 3 , x 4 } and output y ∈ {y 1 , y 2 , y 3 , y 4 } into components S 1 , S 2 , and T for compositional analysis. These components have a probability of being visited of 0.25, 0.4, and 0.35 respectively, accounting for both input distribution and internal randomness. Each component can be visited on specific inputs, and probabilistically produces an output according to the input value. For instance, the component S 2 is chosen with the probability 0.4 and given input x 1 , it outputs y 3 or y 4 with the equal probability 0.5. The division into components can be used to compositionally produce a joint distribution P XY on the input-output behavior of the program.

analyst does not require prior knowledge of the internal structure of the system. The analyst runs the system many times and records the execution traces to construct an approximate model of the system. Even when the formal specification or precise model of the system is not provided to the analyst, statistical analysis can be directly applied to the system if the analyst can execute the black-box implementation. Due to this random sampling of the systems, statistical analysis provides only approximate estimates. However, it can evaluate the precision and accuracy of the analysis for instance by providing the confidence intervals of the estimated values. One of the important challenges in statistical analysis is to estimate entropy-based properties in probabilistic systems. For example, statistical methods [CCG10, CKN13, CKNP13, CKN14, BP14] have been studied for quantitative information flow analysis [CHM01, KB07, Mal07, CPP08], which estimates an entropy-based property to quantify the leakage of confidential information in a system. More specifically, the analysis estimates mutual information or other properties between two random variables on the secrets and on the observable outputs in the system to measure the amount of information that is inferable about the secret by observing the output. The main technical difficulties in the estimation of entropy-based properties are:

1. to efficiently compute large matrices that represent probability distributions, and 2. to provide a statistical method for correcting the bias of the estimate and computing a confidence interval to evaluate the accuracy of the estimation.

To overcome these difficulties, we propose a method for statistically estimating mutual information, one of the most popular entropy-based properties, in discrete systems. The new method, called hybrid statistical estimation method, integrates black-box statistical analysis and white-box precise analysis, exploiting the advantages of both. More specifically, this method employs some prior knowledge on the system and performs precise analysis (e.g., static analysis of the source code or specification) on some components of the system. Since precise analysis computes the exact sub-probability distributions of the components, the hybrid method using precise analysis is more accurate than statistical analysis alone.

Moreover, the new method can combine multiple statistical analyses on different components of the system to improve the accuracy and efficiency of the estimation. This is based on our new theoretical results that extend and generalize previous work [START_REF] Moddemeijer | On estimation of entropy and mutual information of continuous distributions[END_REF][START_REF] Brillinger | Some data analysis using mutual information[END_REF][START_REF] Chatzikokolakis | Statistical measurement of information leakage[END_REF] on purely statistical estimation. As far as we know this is the first work on a hybrid method for estimating entropy-based properties and their confidence intervals. Note that our approach assumes that the system has discrete inputs and outputs and behaves deterministically or probabilistically, while it can also be applied to non-deterministic systems if the non-determinism has been resolved by schedulers.

To illustrate the method we propose, Fig. 1 presents an example of a probabilistic program (Fig. 1a) having input ranging over X = {x 1 , x 2 , x 3 , x 4 } and output over Y = {y 1 , y 2 , y 3 , y 4 }, built up from three overlapping components S 1 , S 2 and T , and the corresponding joint probability distribution P XY (Fig. 1b). To estimate the full joint distribution P XY , the analyst separately computes the joint sub-distribution for the component T by precise analysis, estimates those for S 1 and S 2 by statistical analysis, and then combines these sub-distributions. Since the statistical analysis is based on the random sampling of execution traces, the empirical sub-distributions for S 1 and S 2 are different from the true ones, while the sub-distribution for T is exact. From these approximate and precise sub-distributions, the proposed method can estimate the mutual information for the entire system and evaluate its accuracy by providing a confidence interval. Owing to the combination of different kinds of analyses (with possibly different parameters such as sample sizes), the computation of the bias and confidence interval of the estimate is more complicated than the previous work on statistical analysis.

Contributions

The contributions of this paper are as follows:

• We propose a new method, called hybrid statistical estimation, that combines statistical and precise analysis on the estimation of mutual information (which can also be applied to Shannon entropy and conditional Shannon entropy). Specifically, we show theoretical results on compositionally computing the bias and confidence interval of the estimate from multiple results obtained by statistical and precise analysis; • We present a weighted statistical analysis method with different sample sizes over different components and a method for adaptively optimizing the sample sizes by evaluating the accuracy and cost of the analysis; • We show how to reduce the sample sizes by using prior knowledge about systems, including an abstraction-thensampling technique based on qualitative analysis. In particular, we point out that the state-of-the-art statistical analysis tool LeakWatch [START_REF] Chothia | Leakwatch: Estimating information leakage from java programs[END_REF] incorrectly computes the bias in estimation using the knowledge of the prior distribution, and explain how the proposed approach fixes this problem; • We show that the proposed method can be applied not only to composed systems but also to the source codes of a single system by decomposing it into components and determining the analysis method for each component; • We provide a practical implementation of the method in the HyLeak tool [BKLT], and show how the techniques in this paper can be applied to multiple benchmarks; • We evaluate the quality of the estimation in this method, showing that the estimates are more accurate than statistical analysis alone for the same sample size, and that the new method outperforms the state-of-the-art statistical analysis tool LeakWatch; • We demonstrate the effectiveness of the hybrid method in case studies on the quantification of information leakage.

A preliminary version of this paper, without proofs, appeared in [START_REF] Kawamoto | Hybrid statistical estimation of mutual information for quantifying information flow[END_REF]. Also a preliminary version of the implementation description (Sections 7 and 8.2.1), without details, appeared in the tool paper describing HyLeak [START_REF] Biondi | Hyleak: Hybrid analysis tool for information leakage[END_REF]. In this paper we add the estimation of Shannon entropy (Propositions 4.3, 4.4 and 6.2) and that of conditional entropy (Propositions 5.4 and 5.5). We also show the formulas for the adaptive analysis using knowledge of prior distributions (Proposition 6.3) and using the abstraction-then-sampling technique (Theorem 6.4). Furthermore, we provide detailed explanation on the implementation in the HyLeak tool in Section 7, including how to decompose the source code of a system into components. We also present more experimental results with details in Section 8.2. Finally, we add Appendix A to present the detailed proofs.

The rest of the paper is structured as follows. Section 2 introduces background in quantification of information and compares precise analysis with statistical analysis for the estimation of mutual information. Section 3 overviews our new method for estimating mutual information. Section 4 describes the main results of this paper: the statistical estimation of mutual information for the hybrid method, including the method for optimizing sample sizes for different components. Section 5 presents how to reduce sample sizes by using prior knowledge about systems, including the abstraction-then-sampling technique with qualitative analysis. Section 6 shows the optimal assignment of samples to components to be samples statistically to improve the accuracy of the estimate. Section 7 overviews the implementation of the techniques in the HyLeak tool, including how to decompose the source code of a system into components and to determine the analysis method for each component. Section 8 evaluates the proposed method and illustrates its effectiveness against the state of the art, and Section 9 concludes the paper. Detailed proofs can be found in Appendix A.

Related Work

The information-theoretical approach to program security dates back to the work of Denning [START_REF] Denning | A lattice model of secure information flow[END_REF] and Gray [START_REF] Gray | Toward a mathematical foundation for information flow security[END_REF]. Clark et al. [START_REF] Clark | Quantitative analysis of the leakage of confidential data[END_REF][START_REF] Clark | A static analysis for quantifying information flow in a simple imperative language[END_REF] presented techniques to automatically compute mutual information of an imperative language with loops. For a deterministic program, leakage can be computed from the equivalence relations on the secret induced by the possible outputs, and such relations can be automatically quantified [START_REF] Backes | Automatic discovery and quantification of information leaks[END_REF]. Underand over-approximation of leakage based on the observation of some traces have been studied for deterministic programs [START_REF] Mccamant | Quantitative information flow as network flow capacity[END_REF][START_REF] Newsome | Measuring channel capacity to distinguish undue influence[END_REF]. As an approach without replying on information theory McCamant et al. [START_REF] Kang | DTA++: dynamic taint analysis with targeted controlflow propagation[END_REF] developed tools implementing dynamic quantitative taint analysis techniques for security.

Fremont and Seshia [START_REF] Fremont | Speeding up smt-based quantitative program analysis[END_REF] present a polynomial time algorithm to approximate the weight of traces of deterministic programs with possible application to quantitative information leakage. Progress in randomized program analysis includes a scalable algorithm for uniform generation of sample from a distribution defined as constraints [CFM + 15, CMV13], with applications to constrained-random program verification.

The statistical approach to quantifying information leakage has been studied since the seminal work by Chatzikokolakis et al. [START_REF] Chatzikokolakis | Statistical measurement of information leakage[END_REF]. Chothia et al. have developed this approach in tools leakiEst [START_REF] Chothia | A tool for estimating information leakage[END_REF][START_REF] Chothia | leakiEst[END_REF] and LeakWatch [START_REF] Chothia | Leakwatch: Estimating information leakage from java programs[END_REF][START_REF] Chothia | LeakWatch[END_REF]. The hybrid statistical method in this paper can be considered as their extension with the inclusion of component weighting and adaptive priors inspired by the importance sampling in statistical model checking [START_REF] Barbot | Coupling and importance sampling for statistical model checking[END_REF][START_REF] Clarke | Statistical model checking for cyber-physical systems[END_REF]. To the best of our knowledge, no prior work has applied weighted statistical analysis to the estimation of mutual information or any other leakage measures.

The idea on combining static and randomized approaches to quantitative information flow was first proposed by Köpf and Rybalchenko [START_REF] Köpf | Approximation and randomization for quantitative information-flow analysis[END_REF] while our approach takes a different approach relying on statistical estimation to have better precision and accuracy and is general enough to deal with probabilistic systems under various prior information conditions. In related fields, the hybrid approach combining precise and statistical analysis have been proven to be effective, for instance in concolic analysis [START_REF] Majumdar | Hybrid concolic testing[END_REF][START_REF] Liu | Hybrid test data generation[END_REF], where it is shown that input generated by hybrid techniques leads to greater code coverage than input generated by both fully random and concolic generation. After the publication of the preliminary version [START_REF] Kawamoto | Hybrid statistical estimation of mutual information for quantifying information flow[END_REF] of this paper, a few papers on quantitative information flow combining symbolic and statistical approaches have been published. Malacaria et al. [MKP + 18] present an approach that performs Monte Carlo sampling over symbolic paths while the prior distributions are restricted to be uniform. Sweet et al. [STS + 18] combine abstraction interpretation with sampling and concolic execution for reasoning about Bayes vulnerability. Unlike our work, these two studies aim at giving only bounds on information leakage and do not use statistical hypothesis testing.

Our tool HyLeak processes a simple imperative language that is an extension of the language used in the QUAIL tool version 2.0 [START_REF] Biondi | Comparative analysis of leakage tools on scalable case studies[END_REF]. The algorithms for precise computation of information leakage used in this paper are based on trace analysis [START_REF] Biondi | Quantifying information leakage of randomized protocols[END_REF], implemented in the QUAIL tool [START_REF] Biondi | QUAIL[END_REF][START_REF] Biondi | QUAIL: A quantitative security analyzer for imperative code[END_REF][START_REF] Biondi | Comparative analysis of leakage tools on scalable case studies[END_REF]. As remarked above, the QUAIL tool implements only a precise calculation of leakage that examines all executions of programs. Hence the performance of QUAIL does not scale, especially when the program performs complicated computations that yield a large number of execution traces. The performance of QUAIL as compared to HyLeak is represented by the "precise" analysis approach in Section 8. Since QUAIL does not support the statistical approach or the hybrid approach, it cannot handle large problems that HyLeak can analyze.

As remarked above, the stochastic simulation techniques implemented in HyLeak have also been developed in the tools LeakiEst [START_REF] Chothia | A tool for estimating information leakage[END_REF] (with its extension [START_REF] Kawamoto | Compositionality results for quantitative information flow[END_REF]) and LeakWatch [START_REF] Chothia | Probabilistic point-to-point information leakage[END_REF][START_REF] Chothia | Leakwatch: Estimating information leakage from java programs[END_REF]. The performance of these tools as compared to HyLeak is represented by the "statistical" analysis approach in Section 8.

The tool Moped-QLeak [START_REF] Chadha | Computing information flow using symbolic model-checking[END_REF] computes the precise information leakage of a program by transforming it into an algebraic decision diagram (ADD). As noted in [START_REF] Biondi | Comparative analysis of leakage tools on scalable case studies[END_REF], this technique is efficient when the program under analysis is simple enough to be converted into an ADD, and fails otherwise even when other tools including HyLeak can handle it. In particular, there are simple examples [START_REF] Biondi | Comparative analysis of leakage tools on scalable case studies[END_REF] where Moped-QLeak fails to produce any result but that can be examined by QUAIL and LeakWatch, hence by HyLeak.

Many information leakage analysis tools restricted to deterministic input programs have been released, including TEMU [START_REF] Newsome | Measuring channel capacity to distinguish undue influence[END_REF], squifc [START_REF] Phan | Abstract model counting: a novel approach for quantification of information leaks[END_REF], jpf-qif [START_REF] Phan | Symbolic quantitative information flow[END_REF], QILURA [START_REF] Phan | Quantifying information leaks using reliability analysis[END_REF], nsqflow [VEB + 16], and SHARPPI [START_REF] Weigl | Efficient sat-based pre-image enumeration for quantitative information flow in programs[END_REF]. Some of these tools have been proven to scale to programs of thousands of lines written in common languages like C and Java. Such tools are not able to compute the Shannon leakage for the scenario of adaptive attacks but only compute the min-capacity of a deterministic program for the scenario of one-try guessing attacks, which give only a coarse upper bound on the Shannon leakage. More specifically, they compute the logarithm of the number of possible outputs of the deterministic program, usually by using model counting on a SMT-constraint-based representation of the possible outputs, obtained by analyzing the program. Contrary to these tools, HyLeak can analyze randomized programs 1 and provides a quite precise estimation of the Shannon leakage of the program, not just a coarse upper bound. As far as we know, HyLeak is the most efficient tool that has this greater scope and higher accuracy.

Background

In this section we introduce the basic concepts used in the paper. We first introduce some notions in information theory to quantify the amount of some information in probabilistic systems. Then we compare two previous analysis approaches to quantifying information: precise analysis and statistical analysis.

Quantification of Information

In this section we introduce some background on information theory, which we use to quantify the amount of information in a probabilistic system. Hereafter we write X and Y to denote two random variables, and X and Y to denote the sets of all possible values of X and Y , respectively. We denote the number of elements of a set A by #A. Given a random variable A we denote by E[A] or by A the expected value of A, and by

V[A] the variance of A, i.e., V[A] = E (A -E[A]) 2 .
The logarithms in this paper are to the base 2 and we often abbreviate log 2 to log.

Channels

In information theory, a channel models the input-output relation of a system as a conditional probability distribution of outputs given inputs. This model has also been used to formalize information leakage in a system that processes confidential data: inputs and outputs of a channel are respectively regarded as secrets and observables in the system and the channel represents relationships between the secrets and observables.

A discrete channel is a triple (X , Y, C) where X and Y are two finite sets of discrete input and output values respectively and C is an #X × #Y matrix where each element C[x, y] represents the conditional probability of an output y given an input x; i.e., for each x ∈ X , y∈Y C[x, y] = 1 and 0 ≤ C[x, y] ≤ 1 for all y ∈ Y.

A prior is a probability distribution on input values X . Given a prior P X over X and a channel C from X to Y, the joint probability distribution P XY of X and Y is defined by: P XY [x, y] = P X [x]C[x, y] for each x ∈ X and y ∈ Y.

Shannon Entropy

We recall some information-theoretic measures as follows. Given a prior P X on input X, the prior uncertainty (before observing the system's output Y ) is defined as

H(X) = - x∈X P X [x] log 2 P X [x]
while the posterior uncertainty (after observing the system's output Y ) is defined as

H(X|Y ) = - y∈Y + P Y [y] x∈X P X|Y [x|y] log 2 P X|Y [x|y],
1 Some of these tools, like jpf-qif and nsqflow, present case studies on randomized protocols. However, the randomness of the programs is assumed to have the most leaking behavior. E.g., in the Dining Cryptographers this means assuming all coins produce head with probability 1.

where P Y is the probability distribution on the output Y , Y + is the set of outputs in Y with non-zero probabilities, and P X|Y is the conditional probability distribution of X given Y :

P Y [y] = x ∈X P XY [x , y] P X|Y [x|y] = P XY [x, y] P Y [y] if P Y [y] = 0.
H(X|Y ) is also called the conditional entropy of X given Y .

Mutual Information

The amount of information gained about a random variable X by knowing a random variable Y is defined as the difference between the uncertainty about X before and after observing Y . The mutual information I(X; Y ) between X and Y is one of the most popular measures to quantify the amount of information on X gained by Y :

I(X; Y ) = x∈X ,y∈Y P XY [x, y] log 2 P XY [x, y] P X [x]P Y [y]
where P Y is the marginal probability distribution defined as

P Y [y] = x∈X P XY [x, y].
In the security scenario, information-theoretical measures quantify the amount of secret information leaked against some particular attacker: the mutual information between two random variables X on the secrets and Y on the observables in a system measures the information that is inferable about the secret by knowing the observable. In this scenario mutual information, or Shannon leakage, assumes an attacker that can ask binary questions on the secret's value after observing the system while min-entropy leakage [START_REF] Smith | On the foundations of quantitative information flow[END_REF] considers an attacker that has only one attempt to guess the secret's value.

Mutual information has been employed in many other applications including Bayesian networks [START_REF] Finn | Introduction to Bayesian Networks[END_REF], telecommunications [START_REF] Robert | Information Theory and Reliable Communication[END_REF], pattern recognition [START_REF] Escolano | Information Theory in Computer Vision and Pattern Recognition[END_REF], machine learning [START_REF] David | Information Theory, Inference & Learning Algorithms[END_REF], quantum physics [START_REF] Mark | Quantum Information Theory[END_REF], and biology [START_REF] Adami | Information theory in molecular biology[END_REF]. In this work we focus on mutual information and its application to the above security scenario.

Computing Mutual Information in Probabilistic Systems

In this section we present two previous approaches to computing mutual information in probabilistic systems in the context of quantitative information flow. Then we compare the two approaches to discuss their advantages and disadvantages.

In the rest of the paper a probabilistic system S is defined as a finite set of execution traces such that each trace tr records the values of all variables in S and is associated with a probability P S [tr ]. Note that S does not have non-deterministic transitions. For the sake of generality we do not assume any specific constraints at this moment.

The main computational difficulty in calculating the mutual information I(X; Y ) between input X and output Y lies in the computation of the joint probability distribution P XY of X and Y , especially when the system consists of a large number of execution traces and when the distribution P XY is represented as a large data structure. In previous work this computation has been performed either by the precise approach using program analysis techniques or by the statistical approach using random sampling and statistics.

Precise Analysis

Precise analysis consists of analyzing all the execution traces of a system and determining for each trace tr , the input x, output y, and probability P S [tr ] by concretely or symbolically executing the system. The precise analysis approach in this paper follows the depth-first trace exploration technique presented by Biondi et al. [START_REF] Biondi | Comparative analysis of leakage tools on scalable case studies[END_REF].

To obtain the exact joint probability P XY [x, y] for each x ∈ X and y ∈ Y in a system S, we sum the probabilities of all execution traces of S that have input x and output y, i.e., P XY [x, y] = P S [tr ] tr ∈ S has input x and output y where P S is the probability distribution over the set of all traces in S. This means the computation time depends on the number of traces in the system. If the system has a very large number of traces, it is intractable for the analyst to precisely compute the joint distribution and consequently the mutual information.

In [START_REF] Yasuoka | Quantitative information flow as safety and liveness hyperproperties[END_REF] the calculation of mutual information is shown to be computationally expensive. This computational difficulty comes from the fact that entropy-based properties are hyperproperties [START_REF] Michael | Hyperproperties[END_REF] that are defined using all execution traces of the system and therefore cannot be verified on each single trace. For example, when we investigate the information leakage in a system, it is insufficient to check the leakage separately for each component of the system, because the attacker may derive sensitive information by combining the outputs of different components. More generally, the computation of entropy-based properties (such as the amount of leaked information) is not compositional, in the sense that an entropy-based property of a system is not the (weighted) sum of those of the components.

For this reason, it is inherently difficult to naïvely combine analyses of different components of a system to compute entropy-based properties. In fact, previous studies on the compositional approach in quantitative information flow analysis have faced certain difficulties in obtaining useful bounds on information leakage [BK11, ES13, KG15, KCP17].

Statistical Analysis

Due to the complexity of precise analysis, some previous studies have focused on computing approximate values of entropy-based measures. One of the common approaches is statistical analysis based on Monte Carlo methods, in which approximate values are computed from repeated random sampling and their accuracy is evaluated using statistics. Previous work on quantitative information flow has used statistical analysis to estimate mutual information [START_REF] Chatzikokolakis | Statistical measurement of information leakage[END_REF][START_REF] Moddemeijer | On estimation of entropy and mutual information of continuous distributions[END_REF][START_REF] Brillinger | Some data analysis using mutual information[END_REF], channel capacity [CCG10, BP14] and min-entropy leakage [START_REF] Chothia | Leakwatch: Estimating information leakage from java programs[END_REF][START_REF] Chothia | Statistical estimation of min-entropy leakage[END_REF].

In the statistical estimation of mutual information between two random variables X and Y in a probabilistic system, the analyst executes the system many times and collects the execution traces, each of which has a pair of values (x, y) ∈ X × Y corresponding to the input and output of the trace. This set of execution traces is used to estimate the empirical joint distribution PXY of X and Y and then to estimate the mutual information Î(X; Y ).

Note that the empirical distribution PXY is different from the true distribution P XY and thus the estimated mutual information Î(X; Y ) is different from the true value I(X; Y ). In fact, it is known that entropy-based measures such as mutual information and min-entropy leakage have some bias and variance that depends on the number of collected traces, the matrix size and other factors. However, results on statistics allow us to correct the bias of the estimate and to compute the variance (and the 95% confidence interval). This way we can guarantee the quality of the estimation, which differentiates the statistical approach from the testing approach.

Comparing the Two Analysis Methods

The cost of the statistical analysis is proportional to the size #X × #Y of the joint distribution matrix (strictly speaking, to the number of non-zero elements in the matrix). Therefore, this method is significantly more efficient than precise analysis if the matrix is relatively small and the number of all traces is very large (for instance because the system's internal variables have a large range).

On the other hand, if the matrix is very large, the number of executions needs to be very large to obtain a reliable and small confidence interval. In particular, for a small sample size, statistical analysis does not detect rare events, i.e., traces with a low probability that affect the result. Therefore the precise analysis is significantly more efficient than statistical analysis if the number of all traces is relatively small and the matrix is relatively large (for instance because the system's internal variables have a small range).

The main differences between precise analysis and statistical analysis are summarized in Table 1.

Overview of the Hybrid Statistical Estimation Method

In this section we overview a new method for estimating the mutual information between two random variables X (over the inputs X ) and Y (over the outputs Y) in a system. The method, we call hybrid statistical estimation, integrates both precise and statistical analyses to overcome the limitations on those previous approaches (explained in Section 2.2).

In our hybrid analysis method, we first decompose a given probabilistic system S into mutually disjoint components, which we will define below, and then apply different types of analysis (with possibly different parameters) on different components of the system. More specifically, for each component, our hybrid method chooses the faster analysis between the precise and statistical analyses. Hence the hybrid analysis of the whole system is faster than the precise analysis alone and than the statistical analysis alone, while it gives more accurate estimates than the statistical analysis alone, as shown experimentally in Section 8.

To introduce the notion of components we recall that in Section 2.2 a probabilistic system S is defined as the set of all execution traces such that each trace tr is associated with probability P[tr ]2 . Formally, a decomposition α of S is defined a collection of mutually disjoint non-empty subsets of S: ∅ ∈ α, S = Si∈α S i , and for any S i , S j ∈ α, S i = S j implies S i ∩ S j = ∅. Then each element of α is called a component. In this sense, components are a partition of the execution traces of S. When S is executed, only one of S's components is executed, since the components are mutually disjoint. Hence S i is chosen to be executed with the probability P[S i ].

In decomposing a system we roughly investigate the characteristics of each component's behaviour to choose a faster analysis method for each component. Note that information about a component like its number of traces and the size of its joint sub-distribution matrix can be estimated heuristically before computing the matrix itself. This will be explained in Section 7; before that section this information is assumed to be available. The choice of the analysis method is as follows:

• If a component's behaviour is deterministic, we perform a precise analysis on it.

• If a component's behaviour is described as a joint sub-distribution matrix over small3 subsets of X and Y, then we perform a statistical analysis on the component. • If a component's behaviour is described as a matrix over large 3 subsets of X and Y, then we perform a precise analysis on the component. • By combining the analysis results on all components, we compute the estimated value of mutual information and its variance (and confidence interval). See Section 4 for details. • By incorporating information from qualitative information flow analysis, the analyst may obtain partial knowledge on components and be able to reduce the sample sizes. See Section 5 for details.

See Section 7 for the details on how to decompose a system. One of the main advantages of hybrid statistical estimation is that we guarantee the quality of the outcome by removing its bias and providing its variance (and confidence interval) even though different kinds of analysis with different parameters (such as sample sizes) are combined together.

Another advantage is the compositionality in estimating bias and variance. Since the sampling of execution traces is performed independently for each component, we obtain that the bias and variance of mutual information can be computed in a compositional way, i.e., the bias/variance for the entire system is the sum of those for the components. This compositionality enables us to find optimal sample sizes for the different components that maximize the accuracy of the estimation (i.e., minimize the variance) given a fixed total sample size for the entire system. On the other hand, the computation of mutual information itself is not compositional [START_REF] Kawamoto | On the compositionality of quantitative information flow[END_REF]: it requires calculating the full joint probability distribution of the system by summing the joint sub-distributions of all components of the system.

Finally, note that these results can be applied to the estimation of Shannon entropy (Section 4.3) and conditional Shannon entropy (Section 5.1.3) as special cases. The overview of all results is summarized in Table 2.

Hybrid Method for Statistical Estimation of Mutual Information

In this section we present a method for estimating the mutual information between two random variables X (over the inputs X ) and Y (over the outputs Y) in a system, and for evaluating the precision and accuracy of the estimation. Table 2. Our results on the hybrid method.

We consider a probabilistic system S that consists of (m + k) components S 1 , S 2 , . . . , S m and T 1 , T 2 , . . . , T k each executed with probabilities θ 1 , θ 2 , . . . , θ m and ξ 1 , ξ 2 , . . . , ξ k , i.e., when S is executed, S i is executed with the probability θ i and T j with the probability ξ j . Let I = {1, 2, . . . , m} and J = {1, 2, . . . , k}, one of which can be empty. Then the probabilities of all components sum up to 1, i.e., i∈I θ i + j∈J ξ j = 1. We assume that the analyst is able to compute these probabilities by precise analysis. In the example in Fig. 1 the probabilities of the components S 1 , S 2 , and T are explicitly given as θ 1 = 0.25, θ 2 = 0.45, and ξ 1 = 0.3. However, in general they would be computed by analyzing the behavior of the system before the system executes the three components. More details about how to obtain this in practice are provided when discussing implementation in Section 7.

Once the system is decomposed into components, each component is analyzed either by precise analysis or by statistical analysis. We assume that the analyst can run the component S i for each i ∈ I to record a certain number of S i 's execution traces, and precisely analyze the components T j for j ∈ J to record a certain symbolic representation of T j 's all execution traces, e.g., by static analysis of the source code (or of a specification that the code is known to satisfy). In the example in Fig. 1, this means that the components S 1 and S 2 will be analyzed statistically producing an approximation of their joint distributions, while the component T will be analyzed precisely obtaining its exact joint distribution. The two estimates and one precise joint distributions will be composed to obtain a joint distribution estimate for the whole system, as illustrated in Fig. 1b.

In the rest of this section we present a method for computing the joint probability distribution PXY (Section 3), for estimating the mutual information Î(X; Y ) (Section 4.1), and for evaluating the accuracy of the estimation (Section 4.2). Then we show the application of our hybrid method to Shannon entropy estimation (Section 4.3).

In the estimation of mutual information between the two random variables X and Y in the system S, we need to estimate the joint probability distribution P XY of X and Y .

In our approach this is obtained by combining the joint sub-probability distributions of X and Y for all the components S i 's and T j 's. More specifically, let R i and Q j be the joint sub-distributions of X and Y for the components S i 's and T j 's respectively. Then the joint (full) distribution P XY for the whole system S is defined by:

P XY [x, y] def = i ∈I R i [x, y] + j ∈J Q j [x, y]
for x ∈ X and y ∈ Y. Note that for each i ∈ I and j ∈ J , the sums of all probabilities in the sub-distribution R i and in Q j respectively equal the probabilities θ i (of executing S i ) and ξ j (of executing T j ).

To estimate the joint distribution P XY the analyst computes

• for each j ∈ J , the exact sub-distribution Q j for the component T j by precise analysis on T j , and

• for each i ∈ I, the empirical sub-distribution Ri for S i from a set of traces obtained by executing S i a certain number n i of times.

More specifically, the empirical sub-distribution Ri is constructed as follows. When the component S i is executed n i times, let K ixy be the number of traces that have input x ∈ X and output y ∈ Y. Then n i = x∈X ,y∈Y K ixy . From these numbers K ixy of traces we compute the empirical joint (full) distribution Di of X and Y by:

Di [x, y] def = K ixy n i .
Since S i is executed with probability θ i , the sub-distribution Ri is given by Ri

[x, y] def = θ i Di [x, y] = θiKixy ni .
Then the analyst sums up these sub-distributions to obtain the joint distribution PXY for the whole system S:

PXY [x, y] def = i∈I Ri [x, y] + j∈J Q j [x, y] = i∈I θ i K ixy n i + j∈J Q j [x, y].
Note that R i and Q j may have different matrix sizes and cover different parts of the joint distribution matrix PXY , so they may have to be appropriately padded with zeroes for the summation.

Estimation of Mutual Information and Correction of its Bias

In this section we present our new method for estimating mutual information and for correcting its bias. For each component S i let D i be the joint (full) distribution of X and Y obtained by normalizing R i :

D i [x, y] = Ri[x,y] θi . Let D Xi [x] = y∈Y D i [x, y], D Yi [y] = x∈X D i [x, y], and D = {(x, y) ∈ X × Y : P XY [x, y] = 0}.
Using the estimated joint distribution PXY we can compute the mutual information estimate Î(X; Y ). Note that the mutual information for the whole system is smaller than (or equals) the weighted sum of those for the components, because of its convexity w.r.t. the channel matrix. Therefore it cannot be computed compositionally from those of the components, i.e., it is necessary to compute the joint distribution matrix PXY for the whole system.

Since Î(X; Y ) is obtained from a limited number of traces, it has bias, i.e., its expected value E Î(X; Y ) is different from the true value I(X; Y ). The bias E Î(X; Y ) -I(X; Y ) in the estimation is quantified as follows.

Theorem 4.1 (Mean of estimated mutual information). The expected value E Î(X; Y ) of the estimated mutual information is given by:

E Î(X; Y ) = I(X; Y ) + i∈I θ 2 i 2n i (x,y)∈D ϕ ixy - x∈X + ϕ ix - y∈Y + ϕ iy + O(n -2 i )
where

ϕ ixy = Di [x,y]-Di [x,y] 2 P XY [x,y] , ϕ ix = DXi [x]-DXi [x] 2 P X [x]
and

ϕ iy = DYi [y]-DYi [y] 2 P Y [y]
.

Proof sketch. Here we present only the basic idea. Appendices A.1 and A.2 present a proof of this theorem by showing a more general claim, i.e., Theorem 5.6 in Section 5.2. By properties of mutual information and Shannon entropy, we have:

E Î(X; Y ) -I(X; Y ) = E Ĥ(X) + Ĥ(Y ) -Ĥ(X, Y ) -H(X) + H(Y ) -H(X, Y ) = E Ĥ(X) -H(X) + E Ĥ(Y ) -H(Y ) -E Ĥ(X, Y ) -H(X, Y ) .
Hence it is sufficient to calculate the bias in Ĥ(X), Ĥ(Y ), and Ĥ(X, Y ), respectively. We calculate the bias in Ĥ(X, Y ) as follows. Let f xy (K 1xy , K 2xy , . . . , K mxy ) be the m-ary function defined by:

f xy (K 1xy , K 2xy , . . . , K mxy ) = i∈I θ i K ixy n i + j∈J Q j [x, y] log i∈I θ i K ixy n i + j∈J Q j [x, y] , which equals PXY [x, y] log PXY [x, y]. Let K xy = (K 1xy , K 2xy , . . . , K mxy ).
Then the empirical joint entropy is:

Ĥ(X, Y ) = - (x,y)∈D PXY [x, y] log PXY [x, y] = - (x,y)∈D f xy (K xy ). Let K ixy = E[K ixy ] for each i ∈ I and K xy = E[K xy ].
By the Taylor expansion of f xy (K xy ) (w.r.t. the multiple dependent variables K xy ) at K xy , we have:

f xy (K xy ) = f xy (K xy ) + i∈I ∂fxy (Kxy) ∂Kixy (K ixy -K ixy ) + 1 2 i,j∈I ∂ 2 fxy (Kxy) ∂Kixy∂Kjxy (K ixy -K ixy )(K jxy -K jxy )+ i∈I O(K 3 ixy ).
We use the following properties:

• E K ixy -K ixy = 0, which is immediate from K ixy = E[K ixy ].
• E (K ixy -K ixy )(K jxy -K jxy ) = 0 if i = j, because K ixy and K jxy are independent.

• E (K ixy -K ixy ) 2 = V[K ixy ] = n i D i [x, y](1 -D i [x, y]). Then E Ĥ(X, Y ) = - (x,y)∈D E[ f xy (K xy ) ] = - (x,y)∈D f xy (K xy ) + 1 2 i∈I ∂ 2 fxy (Kxy) ∂Kixy∂Kixy E (K ixy -K ixy ) 2 + O(K 3 ixy ) = - (x,y)∈D f xy (K xy ) + 1 2 i∈I θ 2 i n 2 i P XY [x,y] n i D i [x, y](1 -D i [x, y]) + O(n -2 i ) = H(X, Y ) - i∈I θ 2 i 2n i (x,y)∈D ϕ ixy + O(n -2 i ),
where the derivation of the equalities is detailed in Appendix A. Hence the bias in estimating H(X, Y ) is given by:

E Ĥ(X, Y ) -H(X, Y ) = - i∈I θ 2 i 2n i (x,y)∈D ϕ ixy + O(n -2 i ).
Analogously, we can calculate the bias in Ĥ(X) and Ĥ(Y ) to derive the theorem. See Appendices A.1 and A.2 for the details.

Since the higher-order terms in the formula are negligible when the sample sizes n i are large enough, we use the following as the point estimate of the mutual information:

pe = Î(X; Y ) - i∈I θ 2 i 2n i (x,y)∈D ϕ ixy - x∈X + ϕ ix - y∈Y + ϕ iy
where ϕ ixy , ϕ ix and ϕ iy are respectively empirical values of ϕ ixy , ϕ ix and ϕ iy that are computed from traces; i.e.,

ϕ ixy = Di [x,y]-Di [x,y] 2 PXY [x,y] , ϕ ix = Di [x]-Di [x] 2 PXY [x]
, and

ϕ iy = Di [y]-Di [y] 2 PXY [y]
. Then the bias is closer to 0 when the sample sizes n i are larger.

Evaluation of the Accuracy of Estimation

In this section we present how to evaluate the accuracy of mutual information estimation. The quality of the estimation depends on the sample sizes n i and other factors, and can be evaluated using the variance of the estimate Î(X; Y ). Theorem 4.2 (Variance of estimated mutual information). The variance V Î(X; Y ) of the estimated mutual information is given by:

V Î(X; Y ) = i∈I θ 2 i n i (x,y)∈D D i [x, y] 1 + log P X [x]P Y [y] P XY [x,y] 2 - (x,y)∈D D i [x, y] 1 + log P X [x]P Y [y] P XY [x,y] 2 + O(n -2 i ).
Proof sketch. The variance is calculated using the following:

V Î(X; Y ) = V Ĥ(X) + Ĥ(Y ) -Ĥ(X, Y ) = V Ĥ(X) + V Ĥ(Y ) + V Ĥ(X, Y ) + 2Cov Ĥ(X), Ĥ(Y ) -2Cov Ĥ(X), Ĥ(X, Y ) -2Cov Ĥ(Y ), Ĥ(X, Y ) .
The calculation of these variances and covariances and the whole proof are shown in Appendices A.3 and A.4. (We will present a proof of this theorem by showing a more general claim, i.e., Theorem 5.7 in Section 5.2).

The confidence interval of the estimate of mutual information is also useful to show how accurate the estimate is. A smaller confidence interval corresponds to a more reliable estimate. To compute the confidence interval approximately, we assume the sampling distribution of the estimate Î(X; Y ) as a normal distribution. 4 Then the confidence interval is calculated using the variance v obtained by Theorem 4.2 as follows. Given a significance level α, we denote by z α/2 the z-score for the 100(1 -α 2 ) percentile point. Then the (1 -α) confidence interval of the estimate is given by:

[ max(0, pe -z α/2 √ v), pe + z α/2 √ v ] .
For example, we use the z-score z 0.0025 = 1.96 to compute the 95% confidence interval. To ignore the higher order terms the sample size i∈I n i needs to be at least 4•#X •#Y. By Theorems 4.1 and 4.2, the bias and variance for the whole system can be computed compositionally from those for the components, unlike the mutual information itself. This allows us to adaptively optimize the sample sizes for the components as we will see in Section 6.

Application to Estimation of Shannon Entropy

Hybrid statistical estimation can also be used to estimate the Shannon entropy H(X) of a random variable X in a probabilistic system. Although the results for Shannon entropy are straightforward from those for mutual information, we present the formulas here for completeness. For each i ∈ I let D Xi be the sub-distribution of X for the component S i . Then the mean and variance of the estimate are obtained in the same way as in the Sections 4.1 and 4.2.

Proposition 4.3 (Mean of estimated Shannon entropy). The expected value E Ĥ(X) of the estimated Shannon entropy is given by:

E Ĥ(X) = H(X) - i∈I θ 2 i 2n i x∈X + D Xi [x] (1 -D Xi [x]) P X [x] + O(n -2 i ).
See Appendix A.2 for the proof. From this we obtain the bias of the Shannon entropy estimates.

Proposition 4.4 (Variance of estimated Shannon entropy). The variance V Ĥ(X) of the estimated Shannon entropy is given by:

V Ĥ(X) = i∈I θ 2 i n i x∈X + D Xi [x] 1 + log P X [x] 2 - x∈X + D Xi [x] 1 + log P X [x] 2 + O(n -2 i ).
See Appendix A.4 for the proof. From this we obtain the confidence interval of the Shannon entropy estimates.

Estimation Using Prior Knowledge about Systems

In this section we show how to use prior knowledge about systems to improve the accuracy of the estimation, i.e., to make the variance (and the confidence interval size) smaller and reduce the required sample sizes.

Approximate Estimation Using Knowledge of Prior Distributions

Our hybrid statistical estimation method integrates both precise and statistical analysis, and it can be seen as a generalization and extension of previous work [START_REF] Chatzikokolakis | Statistical measurement of information leakage[END_REF][START_REF] Moddemeijer | On estimation of entropy and mutual information of continuous distributions[END_REF][START_REF] Brillinger | Some data analysis using mutual information[END_REF].

Due to an incorrect computation of the bias, the state-of-the-art statistical analysis tool LeakWatch [CKN14, CKNb] does not correctly estimate mutual information. We explain this problem in Section 5.1.1 and show how to fix it in Section 5.1.2. We extend this result to the estimation of conditional entropy in Section 5.1.3.

State of the Art

For example, Chatzikokolakis et.al. [START_REF] Chatzikokolakis | Statistical measurement of information leakage[END_REF] present a method for estimating mutual information between two random variables X (over secret input values X ) and Y (over observable output values Y) when the analyst knows the (prior) distribution P X of X. In the estimation they collect execution traces by running a system for each secret value x ∈ X . Thanks to the precise knowledge of P X , they have more precise and accurate estimates than the other previous work [START_REF] Moddemeijer | On estimation of entropy and mutual information of continuous distributions[END_REF][START_REF] Brillinger | Some data analysis using mutual information[END_REF] that also estimates P X from execution traces.

Estimation using the precise knowledge of P X is an instance of our result if a system is decomposed into the component S x for each secret x ∈ X = I. If we assume all joint probabilities are non-zero, the following approximate result in [START_REF] Chatzikokolakis | Statistical measurement of information leakage[END_REF] follows from Theorem 4.1.

Corollary 5.1. The expected value E Î(X; Y ) of the estimated mutual information is given by:

E Î(X; Y ) = I(X; Y ) + (#X -1)(#Y-1) 2n + O(n -2 ),
where #X and #Y denote the numbers of possible secrets and observables respectively.

Using this result the bias E Î(X; Y ) -I(X; Y ) is calculated as (#X -1)(#Y-1) 2n in [START_REF] Chatzikokolakis | Statistical measurement of information leakage[END_REF], which depends only on the size of the joint distribution matrix. However, the bias can be strongly influenced by probability values close or equivalent to zero in the distribution; therefore their approximate results can be correct only when all joint probabilities are non-zero and large enough, which is a strong restriction in practice. We show in Section 8.2.3 that the tool LeakWatch [START_REF] Chothia | Leakwatch: Estimating information leakage from java programs[END_REF] uses Corollary 5.1, and consequently miscalculates bias and gives an estimate far from the true value in the presence of very small probability values.

Our Estimation Using Knowledge of Prior Distributions

To overcome these issues we present more general results in the case in which the analyst knows the prior distribution P X . We assume that a system S is decomposed into the disjoint component S ix for each index i ∈ I and input x ∈ X , and that each S ix is executed with probability θ ix in the system S. Let Θ = {θ ix : i ∈ I, x ∈ X }.

Estimation of Mutual Information

In the estimation of mutual information we separately execute each component S ix multiple times to collect execution traces. Unlike the previous work the analyst may change the number of executions

n i P X [x] to n i λ i [x] where λ i [x]
is an importance prior that the analyst chooses to determine how the sample size n i is allocated for each component S ix . An adaptive way of choosing the importance priors will be presented in Section 6. Let Λ = {λ i : i ∈ I}.

Given the number K ixy of S ix 's traces with output y, we define the conditional distribution D i of output given input:

D i [y|x] def = Kixy niλi[x] . Let M ixy = θ 2 ix λi[x] D i [y|x] (1-D i [y|x]
). Then we can calculate the mean and variance of the mutual information ÎΘ,Λ (X; Y ) using Di , Θ, Λ as follows.

Proposition 5.2 (Mean of mutual information estimated using the knowledge of the prior). The expected value E ÎΘ,Λ (X; Y ) of the estimated mutual information is given by:

E ÎΘ,Λ (X; Y ) = I(X; Y ) + i∈I 1 2n i y∈Y + x∈Dy Mixy P XY [x,y] - x∈Dy Mixy P Y [y] + O(n -2 i ).
Proposition 5.3 (Variance of mutual information estimated using the knowledge of the prior).

The variance V ÎΘ,Λ (X; Y ) of the estimated mutual information is given by:

V ÎΘ,Λ (X; Y ) = i∈I x∈X + θ 2 ix n i λ i [x] y∈Dx D i [y|x] log P Y [y] P XY [x,y] 2 - y∈Dx D i [y|x] log P Y [y] P XY [x,y] 2 + O(n -2 i ).
See Appendix A.6 for the proofs.

Estimation of Conditional Entropy

The new method can also estimate the conditional Shannon entropy H(X|Y ) of a random variable X given a random variable Y in a system. In the context of quantitative security, H(X|Y ) represents the uncertainty of a secret X after observing an output Y of the system. The mean and variance of the conditional entropy are obtained from those of the mutual information in the case where the analyst knows the prior. Proof. By H(X|Y ) = H(X) -I(X; Y ), we obtain V ĤΘ,Λ (X|Y ) = V ÎΘ,Λ (X; Y ) . In this section we extend the hybrid statistical estimation method to consider the case in which the analyst knows that the output of some of the components does not depend on the secret input (for instance by static code analysis). Such prior knowledge may help us abstract components into simpler ones and thus reduce the sample size for the statistical analysis.

Abstraction-Then-Sampling Using Partial Knowledge of Components

We illustrate the basic idea of this "abstraction-thensampling" technique as follows. Let us consider an analyst who knows two pairs (x, y) and (x , y ) of inputs and outputs have the same probability in a component S i :

D i [x, y] = D i [x , y ].
Then, when we construct the empirical distribution Di from a set of traces, we can count the number K i{(x,y),(x ,y )} of traces having either (x, y) or (x , y ), and divide it by two:

K ixy = K ix y = K i{(x,y),(x ,y )} 2
. Then to achieve a certain accuracy, the sample size required for the estimation using the prior knowledge on the equality K ixy = K ix y is smaller than that without using it.

In the following we generalize this idea to deal with similar information that the analyst may possess about the components. Let us consider a (probabilistic) system in which for some components, observing the output provides no information on the input. Assume that the analyst is aware of this by qualitative information analysis (for verifying non-interference). Then such a component S i has a sub-channel matrix where all non-zero rows have an identical conditional distribution of outputs given inputs [START_REF] Cover | Elements of information theory[END_REF]. Consequently, when we estimate the #X i × #Y i matrix of S i it suffices to estimate one of the rows, hence the number of executions is proportional to

#Y i instead of #X i × #Y i .
The abstraction-then-sampling approach can be simply explained by referring to the joint distribution matrix in Fig. 2. Note that each row of the sub-distribution matrix for component S 1 is identical, even though the rows of the joint matrix are not, and assume that the analyst knows this by analyzing the code of the program and finding out that for component S 1 the output is independent from the input. Then the analyst would know that it is unnecessary to execute the component separately for each possible input value in S 1 : it is sufficient to execute the component only for one value of the input, and to apply the results to each row in the sub-distribution matrix for component S 1 . This allows the analyst to obtain more precise results and a smaller variance (and confidence interval) on the estimation given a fixed total sample size n i for the component.

Note that even when some components leak no information, computing the mutual information for the whole system requires constructing the matrix of the system, hence the matrices of all components.

Let I be the set of indexes of components that have channel matrices whose non-zero rows consist of the same conditional distribution. For each i ∈ I , we define π i [x] as the probability of having an input x in the component S i .

To estimate the mutual information for the whole system, we apply the abstraction-then-sampling technique to the components I and the standard sampling technique (shown in Section 4) to the components I \ I .

Then the mean and variance of the mutual information are as follows. The following results show that the bias and confidence interval are narrower than when not using the prior knowledge of components.

Theorem 5.6 (Mean of mutual information estimated using the abstraction-then-sampling). The expected value E ÎI (X; Y ) of the estimated mutual information is given by:

E ÎI (X; Y ) = I(X; Y ) + i∈I\I θ 2 i 2n i (x,y)∈D ϕ ixy - x∈X + ϕ ix - y∈Y + ϕ iy + i∈I θ 2 i 2n i (x,y)∈D ψ ixy - y∈Y + ϕ iy + O(n -2 i )
where

ψ ixy def = Di [x,y]πi[x]-Di [x,y] 2 P XY [x,y]
.

See Appendix A.1 for the proof.

Theorem 5.7 (Variance of mutual information estimated using the abstraction-then-sampling). The variance V ÎI (X; Y ) of the estimated mutual information is given by:

V ÎI (X; Y ) = i∈I\I θ 2 i n i (x,y)∈D D i [x, y] 1 + log P X [x]P Y [y] P XY [x,y] 2 - (x,y)∈D D i [x, y] 1 + log P X [x]P Y [y] P XY [x,y] 2 + i∈I θ 2 i n i y∈Y + D Yi [y]γ 2 ixy - y∈Y + D Yi [y]γ ixy 2 + O(n -2 i )
where

γ ixy def = log P Y [y] - x∈X π i [x] log P XY [x, y].
See Appendix A.3 for the proof.

Adaptive Optimization of Sample Sizes

In this section we present a method for deciding the sample size n i of each component S i to estimate mutual information with an optimal accuracy when using the hybrid estimation technique in Section 4 and its variants using prior information on the system in Section 5. The proof for all results in this section can be found in Appendix A.5.

Mutual Information

To decide the sample sizes we take into account the trade-off between accuracy and cost of the statistical analysis: The computational cost increases proportionally to the sample size n i (i.e., the number of S i 's execution traces), while a larger sample size n i provides a smaller variance hence a more accurate estimate. More specifically, given a budget of a total sample size n for the whole system, we obtain an optimal accuracy of the estimate by adjusting each component's sample size n i5 (under the constraint n = i∈I n i ). To compute the optimal sample sizes, we first run each component to collect a small number (compared to n, for instance dozens) of execution traces. Then we calculate certain intermediate values in computing the variance and determine sample sizes for further executions. Formally, let v i be the following intermediate value of the variance for the component S i :

v i = θ 2 i (x,y)∈D Di [x, y] 1 + log P X [x] P Y [y] P XY [x,y] 2 - (x,y)∈D Di [x, y] 1 + log P X [x] P Y [y] P XY [x,y] 2 .
Then we find n i 's that minimize the variance v = i∈I vi ni of the estimate by using the following theorem. Theorem 6.1 (Optimal sample sizes). Given the total sample size n and the above intermediate variance v i of the component S i for each i ∈ I, the variance of the mutual information estimate is minimized if, for all i ∈ I, the sample size n i for S i is given by:

n i = √ vin m j=1 √ vj .
Shannon Entropy Analogously to Theorem 6.1 we can adaptively optimize the sample sizes in the estimation of Shannon entropy in Section 4.3. To compute the optimal sample sizes we define v i by:

v i = θ 2 i x∈X + D Xi [x] 1 + log P X [x] 2 - x∈X + D Xi [x] 1 + log P X [x] 2 + O(n -2 i ).
Then we can compute the optimal sample sizes by using the following proposition.

Proposition 6.2 (Optimal sample sizes for Shannon entropy estimation). Given the total sample size n and the above intermediate variance v i of the component S i for each i ∈ I, the variance of the Shannon entropy estimate is minimized if, for all i ∈ I, the sample size n i for S i satisfies

n i = √ v i n m j=1 √ v j .
Knowledge of the Prior Analogously to Theorem 6.1, the sample sizes n i and the importance priors λ i can be adaptively optimized when the prior distribution of the input is known to the analyst (as presented in Section 5.1).

Proposition 6.3 (Optimal sample sizes when knowing the prior). For each i ∈ I and x ∈ X , let v ix be the following intermediate variance of the component S ix .

v ix = θ 2 ix y∈Dx Di [y|x] log P Y [y] P XY [x,y] 2 - y∈Dx Di [y|x] log P Y [y] P XY [x,y] 2 .
Given the total sample size n, the variance of the estimated mutual information is minimized if, for all i ∈ I and x ∈ X , the sample size n i and the importance prior λ i satisfy:

n i λ i [x] = √ vixn m j=1
√ vjx .

Abstraction-then-sampling Finally, the sample sizes can be optimized for the abstraction-then-sampling approach in Section 5.2 by using the following theorem.

Theorem 6.4 (Optimal sample sizes using the abstraction-then-sampling). Let v i be the following intermediate variance of the component S i :

v i =            θ 2 i (x,y)∈D Di [x, y] 1 + log PX [x] PY [y] PXY [x,y] 2 - (x,y)∈D Di [x, y] 1 + log PX [x] PY [y] PXY [x,y] 2 if i ∈ I \ I θ 2 i y∈Y + D Yi [y]γ 2 ixy - y∈Y + D Yi [y]γ ixy 2 if i ∈ I
Given the total sample size n, the variance of the estimated mutual information is minimized if, for all i ∈ I and x ∈ X , the sample size n i is given by:

n i = √ v i n m j=1 √ v j .

Implementation in the HyLeak Tool

We describe how HyLeak estimates the Shannon leakage of a given program, i.e., the mutual information between secret and output, implementing the hybrid statistical estimation procedure described above. The tool determines which components of the program to analyze with precise analysis and which with statistical analysis, and inserts appropriate annotations in the code. The components are analyzed with the chosen technique and the results are composed into a joint probability distribution of the secret and observable variables. Finally, the mutual information and its confidence interval are computed from the joint distribution.

The HyLeak tool, including user documentation and source code is freely available at https://project.inria.fr/hyleak/. Multiple examples and the scripts to generate the results are also provided.

HyLeak is very simple to use. The code of the system to analyze is written in a file e.g., system.hyleak. We invoke the tool with the command: ./hyleak system.hyleak

The tool generates various .pp text files with analysis information and the control flow graph of the program. Finally, it outputs the prior and posterior Shannon entropy estimates for the secret, the estimated Shannon leakage of the program before and after bias correction, and its confidence interval. HyLeak can also print the channel matrix and additional information; the full list of arguments is printed by ./hyleak -h.

Illustrating Example: Random Walk

Consider the following random walk problem (modeled in Fig. 3).

The secret is the initial location of an agent, encoded by a single natural number representing an approximate distance from a given point, e.g., in meters. Then the agent takes a fixed number of steps. At each step the distance of the agent increases or decreases by 10 meters with the same probability. After this fixed number of random walk steps, the final location of the agent is revealed, and the attacker uses it to guess the initial location of the agent.

This problem is too complicated to analyze by precise analysis, because the analysis needs to explore every possible combination of random paths, amounting to an exponential number in the random walk steps. It is also intractable to analyze with a fully statistical approach, since there are hundreds of possible secret values and the program has to be simulated many times for each of them to sufficiently observe the agent's behavior.

As shown in Section 8, HyLeak's hybrid approach computes the leakage significantly faster than the fully precise analysis and more accurately than the fully statistical analysis.

Architecture

The HyLeak tool implementation consists of the following 4 steps. Steps 1 and 2 are implemented with different ANTLR parsers [START_REF] Parr | The Definitive ANTLR Reference: Building Domain Specific Languages[END_REF]. The implementation of Step 3 inherits some code from the QUAIL tool [START_REF] Biondi | QUAIL[END_REF][START_REF] Biondi | QUAIL: A quantitative security analyzer for imperative code[END_REF][START_REF] Biondi | Comparative analysis of leakage tools on scalable case studies[END_REF] to employ QUAIL's optimization techniques for precise analysis, i.e., parallel analysis of execution traces and compact Markovian state representation.

Step 1: Preprocessing

Step 1a. Lexing, parsing and syntax checking. HyLeak starts by lexical analysis, macro substitution and syntax analysis. In macro substitution the constants defined in the input program are replaced with their declared values, and simple operations are resolved immediately. In the example in Fig. 3, this replaces the value of constant MAX on Line 24 with its declared value from Line 1. The tool checks whether the input program correctly satisfies the language syntax. In case of syntax errors, an error message describing the problem is produced and execution is terminated.

Step 1b. Loop unrolling and array expansion. for loops ranging over fixed intervals are unrolled to optimize the computation of variable ranges and thus program decomposition in Step 2. In the example in Fig. 3, the for loop in Line 24 gets replaced by a fixed number of repetitions of its code with increasing values of the variable time. Similarly, arrays are replaced with multiple variables indexed by their position number in the array. Note that these techniques are used only to optimize program decomposition and not required to compute the leakage in programs with arbitrary loops.

Step 2: Program Decomposition and Internal Code Generation

If a simulate or simulate-abs statement is present in the code, the program decomposition step is skipped and such statements are used to determine program decomposition.

The code may be decomposed heuristically only at conditional branching to efficiently compute the probability of each component being executed. Moreover, each component must be a terminal in the control flow graph, hence no component is executed afterwards. This is because the estimation method requires that the channel matrix for the system is the weighted sum of those for its components, and that the weight of a component is the probability of executing it, as explained in Sections 3. The analysis method and its parameters for each component S i are decided by estimating the computational cost of analyzing S i . Let Z i be the set of all internal randomness (i.e., the non-secret variables whose values are assigned according to probability distributions) in S i . Then the cost of the statistical analysis is proportional to S i 's sub-channel matrix size #X i × #Y i , while the cost of the precise analysis is proportional to the number of all traces in S i (in the worst case proportional to #X i × #Z i ). Hence the cost estimation is reduced to counting #Y i and #Z i .

To obtain this, for each variable and each code line, an estimation of the number of possible values of the variable at the specific code line is computed. This is used to evaluate at each point in the input program whether it would be more expensive to use precise or statistical analysis. These estimations are shown as comments for different lines of the source code in Fig. 3. To reduce the computational cost of the estimation of variable ranges, we apply ad-hoc heuristics to obtain approximate estimates.

After determining the decomposition of the program, HyLeak automatically adds simulate and/or simulate-abs statements in the code to signal which parts of the input program should be analyzed with standard random sampling (Section 4) and with abstraction-then-sampling (Section 5.2) respectively. For instance, since no annotations originally exist in the example source code in Fig. 3, HyLeak adds the simulate-abs statement (written in red) on Line 23. The procedure for decomposition is shown in Fig. 5 and is illustrated in Section 7.3 using the Random Walk example of Fig. 3. While the decomposition procedure is automated, it is a heuristic that does not guarantee to produce an optimal decomposition. Hence for usability the choice of analysis can be controlled by user's annotations on the code.

At the end, the input program is translated into a simplified internal language. Conditional statements and loops (if, for, and while) are rewritten into if-goto statements.

Step 3: Program Analysis In this step the tool analyzes the executions of the program using the two approaches.

Step 3a. Precise analysis. The tool performs a depth-first symbolic execution of all possible execution traces of the input program, until it finds a return, simulate, or simulate-abs statement. When reaching a return statement the tool recognizes the symbolic path as terminated and stores its secret and output values. In the cases of simulate and simulate-abs statements it halts the symbolic path, saves the resulting program state, and schedules it for standard random sampling (Section 4) or for abstraction-then-sampling (Section 5.2), respectively, starting from the saved program state. In the example in Fig. 3, the tool analyzes the code precisely and generates one symbolic path for each of the possible if-elseif-else-end statements from Line 7 to Line 22, so 7 symbolic paths in total (as shown in the control flow graph of the code in Fig. 4). Then each of the 7 symbolic paths meets the simulate-abs statement in Line 23, so they get removed from precise analysis and scheduled for abstraction-then-sampling.

Step 3b. Statistical analysis. The tool performs all the statistical analyses and abstraction-thensampling analyses, using the saved program states from Step 3a as starting point of each component to analyze statistically. The sample size for each simulation is automatically decided by using heuristics to have better accuracy with less sample size, as explained in Sections 4 and 5. The results of each analysis is stored as an appropriate joint probability sub-distribution between secret and observable values. In the example in Fig. 3, each of the 7 symbolic paths scheduled for abstraction-then-sampling gets analyzed with the technique. For each of the symbolic paths HyLeak choses a value of the secret and samples only that one, then applies the results for all secret values of the component. HyLeak recomputes the assignment of samples to the components for each 10% of the total samples, following the optimal sample sizes computed for abstraction-thensampling components in Theorem 6.4.

Step 4: Leakage Estimation

In this step the tool aggregates all the data collected by the precise and statistical analyses (performed in Steps 3) and estimates the Shannon leakage of the input program, together with evaluation of the estimation. This is explained in detail together with the program decomposition in Section 7.3.

1. Build the control flow graph of the system. 2. Mark all possible components based on each conditional branching. Each possible component must be a terminal as explained in Section 4. 3. For each possible component S i , check whether it is deterministic or not (by syntactically checking an occurrence of a probabilistic assignment or a probabilistic function call). If it is, mark the component for precise analysis, since deterministic systems necessarily have a smaller number of traces. 4. For each possible component S i , check whether S i 's output variables are independent of its input variables inside S i (by qualitative information flow). If so, mark that the abstraction-then-sampling technique in Section 5.2 is to be used on the component, meaning that component S i will be sampled on a single secret input value and the results will be applied to all secret values. 5. For each S i , estimate an approximate range size #Z i of its internal variables and #Y i of its observable variables. 6. Looking from the leaves to the root of the control flow graph, estimate the cost of statistical and precise analyses, decide the decomposition into components, and mark each component for the cheaper analysis between the two. (For example, use a heuristics that marks each component S i for precise analysis if #Z i ≤ #X i and for statistical analysis otherwise.) 7. Join together adjacent components if they are marked for precise analysis, or if they are marked for statistical analysis and have the same input and output ranges. 8. For each component, perform precise analysis or statistical analysis as marked. Fig. 5. Procedure for the decomposition of a system into components given its source code described in Section 3. The actual implementation here uses the control flow graph of the system to guide the procedure.

On the Division into Components of the Random Walk Benchmark

In this section we briefly discuss how the Random Walk example in Fig. 3 can be divided into components using the procedure in Fig. 5. The procedure in Fig. 5 shows in more detail how to implement the procedure described in Section 3. In the implementation, the construction of a the control flow graph of the system is used to guide the division into components.

The control flow graph generated by HyLeak is shown in Fig. 4. Note that HyLeak has added a simulate-abs statement to the code, visible in the control flow graph.

The control flow graph in Fig. 4 helps understanding how HyLeak has implicitly divided the program into components. Note that the simulate-abs node marked in red has 7 in-edges. Each of these edges corresponds to a different symbolic path (i.e., a set of execution traces following the same edge), with different possible values for the secret variable sec and the internal variable loc. HyLeak has determined heuristically that at Line 23 the number of possible values of the observable variables (TOT_OBS = 1) is smaller than the number of possible values of internal variables (TOT_INT = 601), hence statistical simulation will be more efficient than precise analysis on these components. Also, in the code after line 23, HyLeak has determined by syntactically analyzing the code that the values of the observable variables do not depend on the secret, hence each row of the sub-channel matrix for each of these components is identical, much like component S 1 in Fig. 2.

Hence, the abstraction-then-sampling technique of Section 5.2 can be applied, meaning that the behavior of each component will be simulated only for a single value of the secret in the set of possible secret values of the component, and the results will be applied to each row of the channel matrix. Now that the analysis has gathered all the necessary information, HyLeak computes the leakage of the system under analysis. More specifically, it constructs an (approximate) joint posterior distribution of the secret and observable values of the input program from all the collected data produced by Step 3, as explained in Section 3. Then the tool estimates the Shannon leakage value from the joint distribution, including bias correction (see Sections 4 and 5). Finally, a 95% confidence interval for the estimated leakage value is computed to roughly evaluate the quality of the analysis.

In the example in Fig. 3, HyLeak outputs the prior Shannon entropy 8.9658, the posterior Shannon entropy 7.0428, the Shannon leakage (after bias correction) 1.9220, and the confidence interval [1.9214, 1.9226].

Evaluation

We evaluate experimentally the effectiveness of our hybrid method compared to the state of the art. We first discuss the cost and quality of the estimation, then test the hybrid method against fully precise/fully statistical analyses on Shannon leakage benchmarks.

On the Tradeoff between the Cost and Quality of Estimation

In the hybrid statistical estimation, the estimate takes different values probabilistically, because it is computed from a set of traces that are generated by executing a probabilistic system. Fig. 6 shows the sampling distribution of the mutual information estimate of the joint distribution in Fig. 1b in Section 1. The graph shows the frequency (on the y axis) of the mutual information estimates (on the x axis) when performing the estimation 1000 times. In each estimation we perform precise analysis on the component T and statistical analysis on S 1 and S 2 (with a sample size of 5000). The graph is obtained from 1000 samples each of which is generated by combining precise analysis on a component and statistical analysis on 2 components (using 5000 randomly generated traces). As shown in Fig. 6 the estimate after the correction of bias by Theorem 4.1 is closer to the true value. The estimate is roughly between the lower and upper bounds of the 95% confidence interval calculated using Theorem 4.2. The interval size depends on the sample size in statistical analysis as shown in Fig. 8a. In Fig. 8a we illustrated the relationships between the size of the confidence interval and the sample size in the statistical analysis. We used an example with the randomly generated 10 × 10 channel matrix presented in Fig. 7 and the uniform prior. The graph shows the frequency (on the y axis) of the corrected mutual information estimates (on the x axis) that are obtained by estimating the mutual information value 1000, 5000 and 10000 times. When the sample size is k times larger then the confidence interval is

√ k times narrower.
The interval size also depends on the amount of precise analysis as shown in Fig. 8b. If we perform precise analysis on larger components, then the sampling distribution becomes more centered (with shorter tails) and the confidence interval becomes narrower. For instance, in Fig. 8b we illustrated the relationships between the size of the confidence interval and the amount of precise analysis. The graph shows the frequency (on the y axis) of the corrected mutual information estimates (on the x axis) that are obtained by estimating the mutual information value 1000 times when statistical analysis is applied to a 10 × 2, 10 × 5 and 10 × 10 sub-matrix of the full 10 × 10 matrix. Using statistical analysis only on a smaller component (10 × 2 sub-matrix) yields a smaller confidence interval than using it on the whole system (10 × 10 matrix). More generally, if we perform precise analysis on larger components, then we have a smaller confidence interval. This means that the hybrid approach produces better estimates than the state of the art in statistical analysis. Due to the combination with precise analysis, the confidence interval estimated by our approach is smaller than LeakWatch [START_REF] Chothia | Leakwatch: Estimating information leakage from java programs[END_REF] for the same sample size.

Shannon Leakage Benchmarks

We compare the performance of our hybrid method with fully precise/statistical analysis on Shannon leakage benchmarks. Our implementations of precise and statistical analyses are variants of the state-of-the art tools QUAIL [START_REF] Biondi | QUAIL: A quantitative security analyzer for imperative code[END_REF][START_REF] Biondi | QUAIL[END_REF] and LeakWatch [START_REF] Chothia | Leakwatch: Estimating information leakage from java programs[END_REF][START_REF] Chothia | LeakWatch[END_REF] respectively. All experiments are performed on an Intel i7-3720QM 2.6GHz eight-core machine with 8GB of RAM running Fedora 21.

Random Walk

We analyze the random walk example presented in Section 7.1 for different values of number of steps M AX. We plot the computation times and the errors in leakage values computed by the three different methods in the graphs presented in Fig 10 . These graph show again that the execution time of precise analysis grows exponentially to the number of steps M AX, while HyLeak and fully randomized analysis do not require much time even for large values of M AX. In the fully randomized analysis the error is always much larger than when using HyLeak. The reservoir sampling problem [START_REF] Vitter | Random sampling with a reservoir[END_REF] consists of selecting K elements randomly from a pool of N > K elements. We quantify the information flow of the commonly-used Algorithm R [START_REF] Vitter | Random sampling with a reservoir[END_REF], shown in Fig 9, for various values of N and K = N/2. In the algorithm, the first K elements are chosen as the sample, then each other element has a probability to replace one element in the sample. We plot the computation times and the errors in leakage values computed by the three different methods in the graphs presented in Fig 11 . It shows that HyLeak hybrid approach performs faster than the full simulation approach and gives more precise results. Compared to the precise analysis, the hybrid approach run faster when increasing the model's complexity.

Reservoir Sampling

1 const N; // number of elements 2 const K; // selection 3 secret array[N ] of int1 s; 4 observable array[K] of int1 r; 5 public int32 j := 0; 6 for i in [0, K-1] do r[i] := s[i] ; 7 for i in [K, N -1] do 8 j := random(0,i); 9 if j<K then r[j] := s[i]; 10 end

Multiple Lying Cryptographers Protocol

The lying cryptographers protocol is a variant of the dining cryptographer multiparty computation protocol [START_REF] Chaum | The dining cryptographers problem: Unconditional sender and recipient untraceability[END_REF] in which a randomly-chosen cryptographer declares the opposite of what they would normally declare, i.e. they lie if they are not the payer, and do not lie if they are the payer. We consider three simultaneous lying cryptographers implementation in which 8 cryptographers run the protocol on three separate overlapping tables A, B and C with 4 cryptographers each. 1 const N :=4; // number of cryptographers at each table 2 const M :=8; // total number of cryptographers / * these bits represent the coin tosses for the three tables * / 3 public array [N ] of int1 coinA, coinB, coinC; 4 public int32 lies; // this is for the liar / * these bits represent the bits declared by the three cryptographers at each table * / 5 public array [N ] of int1 declA, declB, declC;

/ * these are the outputs at each table * / 6 observable int1 outputA := 0; 7 observable int1 outputB := 0; 8 observable int1 outputC := 0; 9 secret int32 h := [0, M ]; // the secret has M+1 possible values 10 lies := random(1, M ); for c in coinA do c := randombit(0.5); 11 for c in coinB do c:=randombit(0.5); 12 for c in coinC do c:=randombit(0.5); The division into components is executed following the principles in Section 7.3. The hybrid approach divides the protocol into 8 components, one for each possible liar. Then each component is analyzed statistically.

13 for i in [0, N -1] do 14 declA[i]:=coinA[i] xor coinA[(i + 1)%N ]; 15 if h == i + 1 then declA[i]:= ! declA[i]; 16 if lies == i + 1 then declA[i]:= ! declA[i]; 17 outputA := outputA xor declA[i]; 18 declB[i]:=coinB[i] xor coinB[(i + 1)%N ]; 19 if h == i + 3 then declB[i]:= ! declB[i];; 20 if lies==i+3 then declB[i]:= ! declB[i];; 21 outputB := outputB xor declB[i]; 22 declC[i] := coinC[i] xor coinC[(i + 1)%N ]; 23 if h == i + 5 then declC[i]:= ! declC[i]; 24 if lies == i + 5 then declC[i]:= ! declC[i];
The results of the experiment are summarized in Table 3. Note that this example has some zeroes in the probabilities in the channel matrix, which makes the tool LeakWatch [START_REF] Chothia | Leakwatch: Estimating information leakage from java programs[END_REF][START_REF] Chothia | LeakWatch[END_REF] compute an incorrect result due to incorrect bias estimation, as explained in Section 5.1. In fact, the leakage value obtained by LeakWatch using Corollary 5.1 is 0.36245, which is far from the correct value of 0.503 in Table 3. On the other hand, the error of our calculation using Proposition 5.2 is 1.87e-4 even with fully statistical analysis.

Shifting Window

1 const N; // number of elements 2 const W; // window size 3 secret int32 sec := [0,N -1]; 4 observable int32 obs; 5 public int32 minS, sizeS, sizeS, minO, sizeO; 6 minS := random(0,N -W -1); 7 if sec ≥ minS then In the Shifting Window example (Fig. 14) the secret sec can take N possible values, and an interval (called a "window") in the secret domain is randomly selected from 1 to W . If the value of the secret is inside the window, then another window is randomly chosen in a similar way and the program outputs a random value from this second window. Otherwise, the program outputs a random value over the secret domain.

In Fig. 13 we present the results of experiments on the shifting window when increasing the size of the secret domain. The execution time of precise analysis grows proportionally to the secret domain size N while HyLeak and fully randomized analysis do not require much time for a larger N . In the fully randomized analysis the error from the true value grows rapidly while in using HyLeak the error is much smaller. The tool HyLeak can analyze programs that terminate only probabilistically. For instance, the program shown in Fig. 15 has a loop that terminates depending on the randomly generated value of the variable rand. No previous work has presented an automatic measurement of information leakage by probabilistic termination, as precise analysis cannot handle non-terminating programs, which typically causes non-termination of the analysis of the program. On the other hand, the stochastic simulation of this program supported in HyLeak terminates after some number of iterations in practice although it may take long for some program executions to terminate. We analyze this model with the hybrid and full simulation approaches only, as the precise analysis does not terminate. The results are given in Table 3. It shows that also with this problem the hybrid approach performs faster than the full simulation approach.

Smart Grid Privacy

A smart grid is an energy network where users (like households) may consume or produce energy. In Fig. 16 we describe a simple model of a smart grid using the HyLeak language. This example is taken from [START_REF] Biondi | Comparative analysis of leakage tools on scalable case studies[END_REF]. The users periodically negotiate with a central aggregator in charge of balancing the total consumption among several users. In practice each user declares to the aggregator its consumption plan. The aggregator sums up the consumptions of the users and checks if it falls within admitted bounds. If not it answers to the users that the consumption is too low or too high by a certain amount, such that they adapt their demand. This model raises some privacy issues as some attacker can try to guess the consumption of a user, and for instance infer whether or not this particular user is at home.

In Fig. 17 we present the experiment results of this smart grid example for different numbers of users. HyLeak takes less time than both fully precise analysis and fully randomized analysis (as shown in the left figure). Moreover, it is closer to the true value than fully randomized analysis especially when the number of users is larger (as shown in the right figure).

Benchmarks results

In Table 3 we show the results of all the benchmarks using fully precise, fully statistical and hybrid analyses, for a sample size of 50000 executions. Timeout is set at 10 minutes.

The results in Table 3 show the superiority of our hybrid approach compared to the state of the art. The hybrid analysis scales better than the precise analysis, since it does not need to analyze every trace of the system. Compared to fully statistical analysis, our hybrid analysis exploits precise analysis on components of the system where statistical estimation would be more expensive than precise analysis. This allows the hybrid analysis to focus the statistical estimation on components of the system where it converges faster, thus obtaining a smaller confidence interval in a shorter time.

Conclusions and Future Work

We have proposed a hybrid statistical estimation method for estimating mutual information by combining precise and statistical analysis, and for compositionally computing the bias and accuracy of the estimate. This naturally extends to the computation of Shannon entropy and conditional Shannon entropy, generalizing previous approaches on computing mutual information. The method automatically decomposes a system into components and determines which type of analysis is better for each component according to the components' properties.

We have also introduced an algorithm to adaptively find the optimal sample sizes for different components in the statistical analysis to minimize their variance and produce a more accurate estimate given a sample size. Moreover, we have presented how to reduce sample sizes by using prior knowledge about systems, including the abstraction-thensampling technique with qualitative analysis. We have shown how to leverage this information on the system to reduce the computation time and error of the estimation.

1 const N :=9; // the total number of users 2 const S; // the number of users we care about 3 const C:=3; // the possible consumptions level 4 const M :=0; // the consumption level of the attacker 5 const LOW T :=2; // the lower threshold 6 const HIGHT :=9; // the upper threshold / * the observable is the order given by the control system * / 7 observable int32 order; We have provided an implementation in the freely available tool HyLeak. We have shown both theoretical and experimental results to demonstrate that the proposed approach and implementation outperform the state of the art.

Future work includes developing theory and tools that extend our hybrid method to the analysis of other properties and integrate further symbolic abstraction techniques into our estimation method. In particular, we are planning to apply our approach to reasoning about uncertainty propagation in probabilistic programs (e.g., [BGP + 16]) by replacing concentration inequalities with our hybrid statistical estimation techniques. Another possible application of the hybrid analysis is to compute information leakage among adaptive agents in the game-theoretic framework [START_REF] Mário | Information leakage games[END_REF][START_REF] Mário | A Game-Theoretic Approach to Information-Flow Control via Protocol Composition[END_REF], in which each agent probabilistically chooses a strategy that is modeled as a component of a channel. Table 3. Shannon leakage benchmark results using the three different methods (precise, full simulation and hybrid).

The results contain the time (in seconds) taken for the analysis, the value of the leakage (in bits), and the error (in bits) compared to the true result, when the true result has been computed with the precise analysis. The result n/a means either that the experiment cannot be performed on this example, which is the case for the precise analysis of the probabilistic terminating loop, or that the error cannot be computed because the precise analysis was not successful.

A. Proofs of Technical Results

In this section we present the detailed proofs of our results.

Hereafter we denote by Q the joint sub-distribution obtained by summing Q j 's:

Q[x, y] def = j ∈J Q j [x, y] .
We write q xy to denote Q[x, y] for abbreviation. Then q xy is the probability that the execution of the system S yields one of T j 's and has input x and output y.

A.1. Proofs for the Mean Estimation Using the Abstraction-Then-Sampling

In this section we present the proof for Theorem 5.6 in Section 5.2, i.e., the result on mean estimation using the abstraction-then-sampling. To show the theorem we present and prove Propositions A.1, A.2, and A.3 below. First, recall that D is defined as the set of pairs consisting of inputs and outputs that appear with non-zero probabilities in the execution of the whole system S:

D = {(x, y) ∈ X × Y : P XY [x, y] > 0} .
Recall also that I = {1, 2, . . . , m}. Let I = {1, 2, . . . , m } for m ≤ m. Proposition A.1 (Mean of joint entropy estimated using the abstraction-then-sampling).

The expected value E ĤI (X; Y ) of the estimated joint entropy is given by:

E ĤI (X, Y ) = H(X, Y ) - i∈I\I θ 2 i 2n i (x,y)∈D ϕ ixy - i∈I θ 2 i 2n i (x,y)∈D ψ ixy + O(n -2 i )
where

ϕ ixy = Di [x,y]-Di [x,y] 2 P XY [x,y] and ψ ixy def = Di [x,y]πi[x]-Di [x,y] 2 P XY [x,y]
.

Proof. We use notations that we have introduced in the previous proofs. For each i ∈ I, let X i be the set of the elements of X that appear with non-zero probabilities in the component S i .

As explained in Section 5.2 we apply the standard sampling technique (shown in Section 4) to the components S i with i ∈ I \ I , and the abstraction-then-sampling technique to the components S i with i ∈ I . We briefly recall the two techniques below.

Using the standard sampling technique we compute the empirical sub-distribution Ri for S i with i ∈ I \ I as follows. The analyst first runs S i a certain number n i of times to obtain the set of execution traces. Let K ixy be the number of traces that have input x ∈ X and output y ∈ Y. Then n i = x∈X ,y∈Y K ixy . From these numbers K ixy of traces we compute the empirical joint (full) distribution Di of X and Y by:

Di [x, y] def = K ixy n i .
Since S i is executed with probability θ i , the sub-distribution Ri is given by Ri [x, y] def = θ i Di [x, y] = θiKixy ni . On the other hand, we use the abstraction-then-sampling sampling technique to compute the empirical subdistribution Ri for S i with i ∈ I as follows. Recall that for each i ∈ I , π i [x] is the probability of having an input x in the component S i . For each i ∈ I all the non-zero rows of S i 's channel matrix are the same conditional distribution; i.e., for each x, x ∈ X i and y ∈ Y i ,

P XY [x,y] πi[x] = P XY [x ,y] πi[x ]
when π i [x] = 0 and π i [x ] = 0. Therefore it is sufficient to estimate only one of the rows. We execute the component S i with an identical input x ∈ X n i times to record the traces. Let K i•y be the number of traces of the component S i that outputs y. Then we define the empirical joint (full) distribution Di of X and Y as:

Di [x, y] def = π i [x]K i•y n i .
Since S i is executed with probability θ i , the sub-distribution Ri is given by: Ri

[x, y] def = θ i Di [x, y] = θ i π i [x]K i•y n i .
Now the empirical joint probability distribution PXY is computed from the above empirical sub-distributions Ri (obtained either by standard sampling or by abstraction-then-sampling) and the exact sub-distributions Q j (obtained by precise analysis):

PXY [x, y] = q xy + i∈I\I θ i K ixy n i + i∈I θ i π i [x]K i•y n i . (1)
Let Let

K xy = (K 1•y , K 2•y , . . . , K m •y , K (m +1)xy , K (m +2
K ixy = E[K ixy ] for each i ∈ I and K xy = E[K xy ]. Then K ixy = n i D i [x, y] = niRi[x,y] θi , and
K i•y = niDi [x,y] πi[x] = niRi[x,y] θiπi[x]
. By the Taylor expansion of f xy (K xy ) (w.r.t. the multiple dependent variables K xy ) at K xy , we have:

f xy (K xy ) = f xy (K xy ) + i∈I\I ∂fxy (Kxy) ∂Kixy (K ixy -K ixy ) + i∈I ∂fxy (Kxy) ∂Ki•y π i [x](K i•y -K i•y ) + i,j∈I\I 1 2 ∂ 2 fxy (Kxy) ∂Kixy∂Kjxy (K ixy -K ixy )(K jxy -K jxy ) + i∈I\I j∈I 1 2 ∂ 2 fxy (Kxy) ∂Kixy∂Kj•y π j [x](K ixy -K ixy )(K j•y -K j•y ) + i,j∈I 1 2 ∂ 2 fxy (Kxy) ∂Ki•y∂Kj•y π i [x]π j [x](K i•y -K i•y )(K j•y -K j•y ) + O(K 3 xy )
To compute the expected value E ĤI (X, Y ) of the estimated joint entropy, it should be noted that:

• E K ixy -K ixy = 0, which is immediate from K ixy = E[K ixy ]. • E K i•y -K i•y = 0, which is immediate from K i•y = E[K i•y ].
• If i = j then E (K ixy -K ixy )(K jxy -K jxy ) = 0, because K ixy and K jxy are independent.

• If i = j then E (K ixy -K ixy )(K j•y -K j•y ) = 0, because K ixy and K j•y are independent.

• If i = j then E (K i•y -K i•y )(K j•y -K j•y ) = 0, because K i•y and K j•y are independent.
• For each i ∈ I \ I , (K ixy : (x, y) ∈ D) follows the multinomial distribution with the sample size n i and the probabilities D i [x, y] for (x, y) ∈ D, therefore

E (K ixy -K ixy ) 2 = V[K ixy ] = n i D i [x, y](1 -D i [x, y]).
• For each i ∈ I , (K i•y : y ∈ Y + ) follows the multinomial distribution with the sample size n i and the probabilities

Di[x,y] πi[x] for (x, y) ∈ D, therefore E (K i•y -K i•y ) 2 = V[K i•y ] = n i D i [x, y] π i [x] 1 - D i [x, y] π i [x] .
Hence the expected value of f xy (K xy ) is given by:

E[f xy (K xy )] = f xy (K xy ) + i∈I\I 1 2 ∂ 2 fxy (Kxy) ∂K 2 ixy E (K ixy -K ixy ) 2 + i∈I πi[x] 2 2 ∂ 2 fxy (Kxy) ∂K 2 i•y E (K i•y -K i•y ) 2 +O(K 3 xy ).
Therefore the expected value of ĤI (X, Y ) is given by:

E ĤI (X, Y ) = H(X, Y ) - (x,y)∈D i∈I\I 1 2 θ 2 i n 2 i P XY [x,y] n i D i [x, y] (1 -D i [x, y]) + i∈I πi[x] 2 2 θ 2 i n 2 i P XY [x,y] n i Di[x,y] πi[x] 1 -Di[x,y] πi[x] + O(n -2 i ) = H(X, Y ) - i∈I\I θ 2 i 2ni (x,y)∈D Di [x,y](1-Di [x,y]) P XY [x,y] - i∈I θ 2 i 2ni (x,y)∈D Di [x,y](πi[x]-Di [x,y]) P XY [x,y] + i∈I O(n -2 i ).
Proposition A.2 (Mean of marginal output entropy estimated using the abstraction-then-sampling). The expected value E ĤI (Y ) of the empirical output entropy is given by:

E ĤI (Y ) = H(Y ) - i∈I θ 2 i 2n i y∈Y + ϕ iy + O(n -2 i ).
Proof. Recall that D is the set of pairs of inputs and outputs with non-zero probabilities, D x = {y: (x, y) ∈ D} and D y = {x: (x, y) ∈ D}. For each i ∈ I \ I and y ∈ Y let L i•y = x∈Dy K ixy . Recall the empirical joint distribution PXY in Equation (1) in the proof of Proposition A.1. Now the empirical marginal distribution PY on outputs is given by: Let Y + be the set of outputs with non-zero probabilities. Then the empirical marginal entropy is:

PY [y] = x∈Dy PXY [x, y] = x∈Dy q xy + i∈I\I θ i K ixy n i + i∈I θ i π i [x]K i•y n i = x∈Dy q xy + i∈I\I θ i L i•y n i + i∈I θ i K i•y n i . Let K y = (L 1•y , L 2•y , . . . , L m •y , K m +1•y , K m +2•y , . . . , K m•y ),
ĤI (Y ) = - y∈Y + PY [y] log PY [y] = - y∈Y + f y (K y ). Let L i•y = E[L i•y ] for each i ∈ I \ I , and K y = E[K y ]. Then L i•y = x∈Dy K ixy .
By the Taylor expansion of f y (K y ) (w.r.t. the multiple dependent variables K y ) at K y , we have:

f y (K y ) = f y (K y ) + i∈I\I ∂fy (Ky) ∂Li•y (L i•y -L i•y ) + i∈I ∂fy (Ky) ∂Ki•y (K i•y -K i•y ) + i,j∈I\I 1 2 ∂ 2 fy (Ky) ∂Li•y∂Lj•y (L i•y -L i•y )(L j•y -L j•y ) + i∈I\I j∈I 1 2 ∂ 2 fy (Ky) ∂Li•y∂Kj•y (L i•y -L i•y )(K j•y -K j•y ) + i,j∈I 1 2 ∂ 2 fy (Ky) ∂Ki•y∂Kj•y (K i•y -K i•y )(K j•y -K j•y ) + O(K 3 y ) Recall that D Yi [y] = x∈X D i [x, y].
To compute the expected value E ĤI (Y ) of the estimated marginal entropy, it should be noted that:

• E L i•y -L i•y = 0, which is immediate from L i•y = E[L i•y ]. • E K i•y -K i•y = 0, which is immediate from K i•y = E[K i•y ].
• If i = j then E (L i•y -L i•y )(L j•y -L j•y ) = 0, because L i•y and L j•y are independent.

• E (L i•y -L i•y )(K j•y -K j•y ) = 0, because L i•y and K j•y are independent.

• If i = j then E (K i•y -K i•y )(K j•y -K j•y ) = 0, because K i•y and K j•y are independent.

• For i ∈ I \ I , (L i•y : y ∈ Y + ) follows the multinomial distribution with the sample size n i and the probabilities D Yi [y] for y ∈ Y + , therefore

E (L i•y -L i•y ) 2 = V[L i•y ] = n i D Yi [y] (1 -D Yi [y]
) .

• For i ∈ I , (K i•y : y ∈ Y + ) follows the multinomial distribution with the sample size n i and the probabilities D Yi [y] for y ∈ Y + , therefore

E (K i•y -K i•y ) 2 = V[K i•y ] = n i D Yi [y] (1 -D Yi [y]) .
Hence the expected value of f y (K y ) is given by:

E[f y (K y )] = f y (K y ) + i∈I\I 1 2 ∂ 2 fy (Ky) ∂L 2 i•y E (L i•y -L i•y ) 2 + i∈I 1 2 ∂ 2 fy (Ky) ∂K 2 i•y E (K i•y -K i•y ) 2 + O(K 3 y ).
Therefore the expected value of ĤI (Y ) is given by:

E ĤI (Y ) = H(Y ) - y∈Y + i∈I\I θ 2 i 2n 2 i P Y [y] E (L i•y -L i•y ) 2 + i∈I θ 2 i 2n 2 i P Y [y] E (K i•y -K i•y ) 2 + O(n -2 i ) = H(Y ) - i∈I θ 2 i 2n i y∈Y + D Yi [y] (1 -D Yi [y]) P Y [y] + O(n -2 i ).
Proposition A.3 (Mean of marginal input entropy estimated using the abstraction-then-sampling). The expected value E ĤI (X) of the empirical input entropy is given by:

E ĤI (X) = H(X) - i∈I\I θ 2 i 2n i x∈X + ϕ ix + O(n -2 i ).
Proof. For the components i ∈ I , the prior π i [x] is known to the analyst and used in the abstraction-then-sampling technique. Hence these components produce no bias in estimating H(X). For the components i ∈ I \ I , we derive the bias in a similar way to the proof of Proposition A.2. Hence the theorem follows. + O(n -2 i ).

Therefore we obtain the theorem.

A.2. Proofs for the Mean Estimation Using Only the Standard Sampling

In this section we present the proofs for Theorem 4.1 in Section 4.1 and Proposition 4.3 in Section 4.3. .

Proof. If I = ∅, then E Î(X; Y ) = E ÎI (X; Y ) . Hence the claim follows from Theorem 5.6.

Proposition 4.3 (Mean of estimated Shannon entropy). The expected value E Ĥ(X) of the estimated Shannon entropy is given by:

E Ĥ(X) = H(X) - i∈I θ 2 i 2n i x∈X + D Xi [x] (1 -D Xi [x]) P X [x] + O(n -2 i ).
Proof. If I = ∅, then E Ĥ(X) = E ĤI (X) . Hence the claim follows from Proposition A.3.

A.3. Proof for the Variance Estimation Using the Abstraction-Then-Sampling

In this section we present the proof for Theorem 5.7 in Section 5.2, i.e., the result on variance estimation using the abstraction-then-sampling.

To show the proofs we first calculate the covariances between random variables in Lemmas A.4, A.5, A.6, A.7, and A.8 as follows. Recall that the covariance Cov [A, B] between two random variables A and B is defined by: Proof. The claim is immediate from the fact that L i•y and K i •y are independent.

Cov [A, B]
Lemma A.8 (Covariance between K i•y and K i •y ). For any i, i ∈ I , Cov [K i•y , K i •y ] is given by:

Cov [K i•y , K i •y ] =    0 if i = i n i D Yi [y](1 -D Yi [y]) if i = i and y = y -n i D Yi [y]D Yi [y ]
otherwise.

Proof. Let i, i ∈ I . If i = i then K i•y and K i •y are independent, hence their covariance is 0. Otherwise, the claim follows from that fact that for each i ∈ I , (K i•y : y ∈ Y + i ) follows the multinomial distribution with the sample size n i and the probabilities D Yi [y] for y ∈ Y + i . Theorem 5.7 (Variance of mutual information estimated using the abstraction-then-sampling). The variance V ÎI (X; Y ) of the estimated mutual information is given by: V ÎI (X; Y ) = variance of the component S i : Given the total sample size n, the variance of the estimated mutual information is minimized if, for all i ∈ I and x ∈ X , the sample size n i is given by:

v i =            θ 2 i (x,
n i = √ v i n m j=1 √ v j .
Proof. By Proposition A.9 the variance v = i∈I v i ni of estimated mutual information is minimised when

n i = √ v i n m j=1 √ v j
. Hence the theorem follows.

Theorem 6.1 (Optimal sample sizes). Given the total sample size n and the above intermediate variance v i of the component S i for each i ∈ I, the variance of the mutual information estimate is minimized if, for all i ∈ I, the sample size n i for S i is given by: n i = Given the total sample size n, the variance of the estimated mutual information is minimized if, for all i ∈ I and x ∈ X , the sample size n i and the importance prior λ i satisfy:

n i λ i [x] = √ vixn m j=1
√ vjx .

Proof. By Proposition A.9 the variance v = i∈I x∈X + vix ni of estimated mutual information is minimised when

n i λ i [x] = √ vixn m j=1
√ vjx . Hence the theorem follows.

A.6. Proof for the Estimation Using the Knowledge of Priors

In this section we present the proof for Propositions 5.2 and 5.3 in Section 5.1.2, i.e., the result on variance estimation using the knowledge of the prior. + O(n -2 i ).

  A probabilistic program composed of 3 components.

  Joint distribution P XY composed of 3 components.

Proposition 5. 4 (

 4 Mean of estimated conditional entropy). The expected value E ĤΘ,Λ (X|Y ) of the estimated conditional Shannon entropy is given by H(X) -E ÎΘ,Λ (X; Y ) where E ÎΘ,Λ (X; Y ) is the expected value of the mutual information in the case where the analyst knows the prior (shown in Proposition 5.2).Proof. By H(X|Y ) = H(X) -I(X; Y ), we obtain E ĤΘ,Λ (X|Y ) = H(X) -E ÎΘ,Λ (X; Y ) . Proposition 5.5 (Variance of estimated conditional entropy). The variance V ĤΘ,Λ (X|Y ) of the estimated conditional Shannon entropy coincides with the variance V ÎΘ,Λ (X; Y ) of the mutual information in the case where the analyst knows the prior (shown in Proposition 5.3).

Fig. 2 :

 2 Fig.2: Joint distribution composed of 3 components. All the rows of component S 1 are identical, hence abstraction-then-sampling can be used on it.

1

  Fig. 3. Source code for the Random Walk illustrative example explained in Section 7.1. The comments show the estimates for the value ranges of some variables, as computed by HyLeak following Step 2 of Section 7.2, where TOT_OBS represent the estimate of the possible combinations of values of all observable variables and TOT_INT the estimate of the possible combinations of values of all internal variables. The red simulate-abs statement in Line 23 shows where the statement will be automatically added to implement the division into components, as explained in Step 2 of Section 7.2 and in more details in Section 7.3.

Fig. 4 :

 4 Fig. 4: Control flow graph for the input code of Fig. 3. The red node corresponds to the simulate-abs statement on Line 23 of Fig. 3. Each arrow entering the red node comes from a symbolic path with different secret and internal values, inducing a different component.

Fig. 6 .

 6 Fig.6. Distribution of mutual information estimate and its confidence interval.

Fig. 7 .

 7 Fig. 7. Channel matrix for the experiments in Section 8.1.

  Estimates and sample sizes.

  Estimates and the ratio of precise analysis.

Fig. 8 .

 8 Fig.8. Smaller intervals when increasing the sample size or the ratio of precise analysis.

Fig. 9 :

 9 Fig. 9: Reservoir sampling.

Fig. 11 .

 11 Fig. 10. Random walk experimental results.

  Fig. 12. Multiple lying cryptographers.

Fig. 13 .

 13 Fig. 13. Shifting window experimental results.

8

  Fig. 14: Shifting Window.

Fig. 15 :

 15 Fig. 15: Probabilistically Terminating Loop.

Fig. 17 .

 17 Fig. 16. Smart Grid Example.

  )xy , . . . , K mxy ), and f xy (K xy ) be the m-ary function: f xy (K xy ) = q xy + which equals PXY [x, y] log PXY [x, y]. Then the empirical joint entropy is: ĤI (X, Y ) = -(x,y)∈D PXY [x, y] log PXY [x, y] = -(x,y)∈D f xy (K xy ).

,

  and f y (K y ) be the following m-ary function:f y (K y ) =x∈Dy which equals PY [y] log PY [y].

Theorem 5. 6 (= 2 i

 62 Mean of mutual information estimated using the abstraction-then-sampling). The expected value E ÎI (X; Y ) of the estimated mutual information is given by:E ÎI (X; Y ) = I(X; Y ) + Di [x,y]πi[x]-Di [x,y] 2 P XY [x,y] .Proof. By Propositions A.1, A.2, and A.3, we obtain the expected value of the estimated mutual information:E ÎI (X; Y ) = E ĤI (X) + E ĤI (Y ) -E ĤI (X, Y ) = I(X; Y ) + i∈I\I θ 2n i (x,y)∈D Di [x,y]-Di [x,y] 2 P XY [x,y] Di [x,y]πi[x]-Di [x,y] 2 P XY [x,y] -y∈Y + DYi [y]-DYi [y] 2 P Y [y]

Theorem 4. 1 ( 2 i

 12 Mean of estimated mutual information). The expected value E Î(X; Y ) of the estimated mutual information is given by:E Î(X; Y ) = I(X; Y ) + i∈I θ 2n i (x,y)∈D ϕ ixy -x∈X + ϕ ix -y∈Y + ϕ iy + O(n -2 i )whereϕ ixy = Di [x,y]-Di [x,y] 2 P XY [x,y] , ϕ ix = DXi [x]-DXi [x] 2 P X [x]and ϕ iy = DYi [y]-DYi [y] 2 P Y [y]

  def = E[(A -E[A])(B -E[B])] .

Lemma A. 7 (

 7 Covariance between L i•y and K i •y ). For any i ∈ I \ I and any i ∈ I , Cov [L i•y , K i •y ] = 0.

i∈I\I θ 2 iD

 2 n i (x,y)∈D D i [x, y] 1 + log P X [x]P Y [y] P XY [x,y] 2 -(x,y)∈D D i [x, y] 1 + log P X [x]P Y [y] P XY [x,y] + D Yi [y]γ 2 ixy -y∈Y + D Yi [y]γ ixy 2 + O(n -2 i )whereγ ixy def = log P Y [y] -x∈X π i [x] log P XY [x, y].Proof. We first define B ixy and B i•y by the following:• B ixy def = ∂fxy (Kxy) ∂Kixy = θi ni (1 + log P XY [x, y]). • For each i ∈ I \ I , B i•y def = ∂fy (Ky) ∂Li•y = θi ni (1 + log P Y [y]). • For each i ∈ I , B i•y def = ∂fy (Ky) ∂Ki•y = θi ni (1 + log P Y [y]).Then the variance of ĤI (X, Y ) is obtained from Lemmas A.4 and A.8 as follows:V ĤI (X, Y ) = E ĤI (X, Y ) 2 -E ĤI (X, Y ) 2 = i,i ∈I\I (x,y)∈D (x ,y )∈D B ixy B i x y Cov [K ixy , K i x y ] + i,i ∈I (x,y)∈D (x ,y )∈D π i [x]B ixy π i [x ]B i x y Cov [K i•y , K i •y ] + O(n -2 i ) = i∈I\I (x,y)∈D n i B 2 ixy D i [x, y](1 -D i [x, y]) -(x ,y )∈D\{(x,y)} B ixy B ix y D i [x, y]D i [x , y ] + i∈I (x,y)∈D n i π i [x] 2 B 2 ixy D Yi [y](1 -D Yi [y]) -(x ,y )∈D\{(x,y)} π i [x]π i [x ]B ixy B ix y D Yi [y]D Yi [y ] D i [x, y]B ixy B ixy -(x ,y )∈D B ix y D i [x , y ] Yi [y] x∈Dy π i [x] 2 B 2 ixy -(x,y)∈D D Yi [y]π i [x]B ixy 2 + O(n -2 i ).

  y)∈D Di [x, y] 1 + log PX [x] PY[y]

  vj . Proof. Let I = ∅. Then this theorem immediately follows from Theorem 6.4. Proposition 6.2 (Optimal sample sizes for Shannon entropy estimation). Given the total sample size n and the above intermediate variance v i of the component S i for each i ∈ I, the variance of the Shannon entropy estimate is minimized if, for all i ∈ I, the sample size n i for S i satisfies n i = √ By Proposition A.9 the variance v = i∈I v i ni of estimated Shannon entropy is minimised when n i = √ follows. Proposition 6.3 (Optimal sample sizes when knowing the prior). For each i ∈ I and x ∈ X , let v ix be the following intermediate variance of the component S ix . v ix = θ 2 ix y∈Dx Di [y|x] log P Y [y] P XY [x,y] 2 -y∈Dx Di [y|x] log P Y [y] P XY [x,y] 2 .

  Proposition 5.2 (Mean of mutual information estimated using the knowledge of the prior). The expected value E ÎΘ,Λ (X; Y ) of the estimated mutual information is given by: E ÎΘ,Λ (X; Y ) = I(X; Y ) + i∈I 1 2n i y∈Y + x∈Dy Mixy P XY [x,y] -x∈Dy Mixy P Y [y]

Table 1 .

 1 Comparison of the precise and statistical analysis methods.

		Precise analysis	Statistical analysis
	Type	White box	Black/gray box
	Analyzes	Source code	Implementation
	Produces	Exact value	Estimate & accuracy evaluation
	Reduces costs by	Random sampling	Knowledge of code & abstraction
	Impractical for	Large number of traces Large channel matrices

  Table A hosts cryptographers 1 to 4, Table B hosts cryptographers 3 to 6, and Table C hosts cryptographers 5 to 8. The identity of the payer is the same in all tables.

Note that this work considers only probabilistic systems without non-deterministic transitions.

Relatively to the number of all execution traces of the component.

In fact, Brillinge[START_REF] Brillinger | Some data analysis using mutual information[END_REF] shows that the sampling distribution of mutual information values is approximately normal for large sample size, and Chatzikokolakis et al.[START_REF] Chatzikokolakis | Statistical measurement of information leakage[END_REF] employ this fact to approximately compute the confidence interval. We also empirically verified that the sampling distribution is closer to the normal distribution when the n i 's are larger enough, and the evaluation of the obtained confidence interval will be demonstrated by experiments in Section 8.1.

This idea resembles the importance sampling in statistical model checking in that the sample size is adjusted to make the estimate more accurate.

obs := BOU N D;

Note that the abstraction-then-sampling relies on partial knowledge on the prior, i.e., the analyst knows π i [x] for all i ∈ I , hence V ĤI (X) has no term for I . On the other hand, the standard sampling here does not use knowledge on the prior.

Acknowledgment

This work was supported by JSPS KAKENHI Grant Number JP17K12667, by JSPS and Inria under the Japan-France AYAME Program, by the MSR-Inria Joint Research Center, by the Sensation European grant, and by région Bretagne.

Probabilistically Terminating Loop

1 const N ; // number of secrets 2 const BOU N D; 3 secret int32 sec := [0, N ]; 4 observable int32 obs; 5 public int32 time := 0; 6 public int2 terminate := 0; 7 public int32 rand; 8 while terminate = 1 do Lemma A.4 (Covariance between K ixy and K i x y ). For any i, i ∈ I \ I , Cov [K ixy , K i x y ] is given by:

Proof. Let i, i ∈ I \ I . If i = i then K ixy and K i x y are independent, hence their covariance is 0. Otherwise, the theorem follows from that fact that for each i ∈ I \ I , (K ixy : (x, y) ∈ D) follows the multinomial distribution with the sample size n i and the probabilities D i [x, y] for (x, y) ∈ D.

Lemma A.5 (Covariance between K ixy and L i •y ). For any i, i ∈ I \ I , Cov [K ixy , L i •y ] is given by:

Proof. Let i, i ∈ I \ I . If i = i then K ixy and L i •y are independent, hence their covariance is 0. Otherwise, the covariance Cov [K ixy , L i•y ] is calculated as:

x ∈Dy

Hence, when y = y :

When y = y :

Lemma A.6 (Covariance between L i•y and L i •y ). For any i ∈ I \ I , Cov [L i•y , L i•y ] is given by:

Proof. Let i, i ∈ I \ I . If i = i then L i•y and L i •y are independent, hence their covariance is 0. Otherwise, the covariance is calculated as:

x ∈Dy

Hence, when y = y :

When y = y :

The variances of ĤI (Y ) is obtained from Lemmas A.6 and A.8 as follows:

Similarly, for

, the variance of ĤI (X) is given by:

which is symmetric to V ĤI (Y ) only w.r.t. I \ I 6 .

The covariance between ĤI (X, Y ) and ĤI (Y ) is obtained from Lemmas A.5 and A.8 as follows:

Similarly, the covariance between ĤI (X, Y ) and ĤI (X) is given by:

The covariance between ĤI (X) and ĤI (Y ) is given by:

Therefore the variance of the mutual information is as follows:

A.4. Proof for the Variance Estimation Using Only the Standard Sampling

In this section we present the proofs for Theorem 4.2 in Section 4.2 and Proposition 4.4 in Section 4.3.

Theorem 4.2 (Variance of estimated mutual information). The variance V Î(X; Y ) of the estimated mutual information is given by:

Proof.

Hence the claim follows from Theorem 5.7.

Proposition 4.4 (Variance of estimated Shannon entropy). The variance V Ĥ(X) of the estimated Shannon entropy is given by:

Then by the proof of Theorem 5.7 in Appendix A.3, we have:

A.5. Proofs for Adaptive Analysis

In this section we present the proofs for the results in Section 6. To prove these it suffices to show the following proposition:

Proposition A.9. Let v 1 , v 2 , . . . , v m be m positive real numbers. Let n, n 1 , n 2 , . . . , n m be (m + 1) positive real numbers such that

The equality holds when

Proof. The proof is by induction on m. When m = 1 the equality holds trivially. When m = 2 it is sufficient to prove

By n 1 , n 2 > 0, this is equivalent to

. We obtain this by:

Next we prove the inductive step as follows.

Finally, when

√ vj for all i = 1, 2, . . . , m, then

Theorem 6.4 (Optimal sample sizes using the abstraction-then-sampling). Let v i be the following intermediate Proof. Since the precise prior is provided to the analyst, we have

By using the results on E Ĥ(Y ) and E Ĥ(X, Y ) in the proof of Theorem 4.1, we obtain the proposition.

Proposition 5.3 (Variance of mutual information estimated using the knowledge of the prior).

The variance V ÎΘ,Λ (X; Y ) of the estimated mutual information is given by: V ÎΘ,Λ (X; Y ) = Proof. Since the precise prior is provided to the analyst, we have

By using the results on V Ĥ(Y ) , V Ĥ(X, Y ) Cov Ĥ(Y ), Ĥ(X, Y ) in the proof of Theorem 4.2, we obtain the proposition.