Once is Never Enough:
Foundations for Sound Statistical Inference in Tor Network Experimentation

ExperimenTor [8] is a Tor experimentation testbed built on top of the ModelNet [72] network emulation platform. ExperimenTor consists of two components that generally run on independent servers (or clusters): one component runs client processes and the other runs the ModelNet core emulator that connects the processes in a virtual network topology. The performance of this architecture was improved in SNEAC [63] through the use of Linux Containers and the kernel’s network emulation module netem, while tooling and network orchestration were improved in NetMirage [70].

Shadow [29] is a hybrid discrete-event network simulator that runs applications as plugins. We provide more background on Shadow in §4.1. Shadow’s original design was improved with the addition of a user-space non-preemptive thread scheduler [51], and later with a high performance dynamic loader [69]. Additional contributions have been made through several research projects [35, 37, 38], and we make further contributions that improve Shadow’s efficiency and correctness as described in §4.2.

We produce more representative networks by considering the state of the network over time rather than modeling a single snapshot as did previous work [32, 38]. We consider relay churn to demonstrate how the true Tor network changes over time. Figure 1 shows the rate of relay churn over all 744 consensus files (1 per hour) in Tor during January 2019. After 2 weeks, fewer than 75% of relays that were part of the network on 2019-01-01 still remain while more than 3,000 new relays joined the network. After 3 weeks, more new relays had joined the network than had remained since 2019-01-01. Our models account for churn by sampling from all such relays as described in §3.2.

In addition to producing more representative models, our Shadow network stack enhancements further improve network accuracy. To demonstrate these contributions, we simulate ten Tor network models that were generated following the methods in §3.2 (using Tor network state from 2019-01). We model Tor at the same $s=0.31$ scale that was used in previous work [38] (i.e., $\approx$2k relays and $\approx$250k users) using a process scale factor of $p=0.01$ (i.e., each TGen process simulated $1/0.01=100$ Tor users). We compare our simulation results to those produced by state-of-the-art methods [38] (which used Tor network state from 2018-01) and to reproducible Tor metrics [66, 67] from the corresponding modeling years (2019 for our work, 2018 for the CCS 2018 work).¹²¹²12Although the models used Tor data spanning one month, we consider it reasonable to reflect the general state of Tor throughout the respective year.

The results in Figure 2 generally show that previous work is noticeably less accurate when compared to Tor 2018 than our work is compared to Tor 2019. We notice that previous work exhibited a high client download error rate in Figure 2(c) and significantly longer download times in Figures 2(e)–2(g) despite the network being appropriately loaded as shown in Figure 2(h). We attribute these errors to the connection limit and network stack limitations that were present in the CCS 2018 version of Shadow (the errors are not present in this work due to our Shadow improvements from §4.2). Also, we remark that the relay goodput in Figure 2(h) exhibits more variance in Tor than in Shadow because the Tor data is being aggregated over a longer time period (1 year for Tor vs. less than 1 hour for Shadow) during which the Tor network composition is significantly changing (see Figure 1).

Our new models and Shadow enhancements enable researchers to run larger networks faster than was previously possible. We demonstrate our improvements to scalability in two ways. First, we compare in the top part of Table 2 the resources required for the 31% experiments described above. We distinguish total run time from the time required to bootstrap all Tor relays and clients, initialize all traffic generators, and reach steady state. We reduced the time required to execute the bootstrapping process by 2 days, 18 hours, or 80%, while we reduced the total time required to run the bootstrapping process plus 25 simulated minutes of steady state by 33 days, 12 hours, or 94%. The ratio of real time units required to execute each simulated time unit during steady state (i.e., after bootstrapping has completed) was reduced by 96%, further highlighting our achieved speedup. When compared to models of the same $s=31\%$ scale from previous work, we observed that our improvements reduced the maximum RAM required to run bootstrapping plus 25 simulated minutes of steady state from 2.6 TiB down to 932 GiB (a total reduction of 1.7 TiB, or 64%).

Second, we demonstrate how our improvements enable us to run significantly larger models by running three Tor models at scale $s=1.0$, i.e., at 100% of the size of the true Tor network. We are the first to simulate Tor test networks of this scale.¹³¹³13We attempted to run a 100% scale Tor network using the CCS 2018 model [38], but it did not complete the bootstrapping phase within 30 days. The bottom part of Table 2 shows that each of our 100% Tor networks consumed at most 3.9 TiB of RAM, completed bootstrapping in 2 days, 21 hours, and ran the entire simulation (bootstrapping plus 25 simulated minutes of steady state) in 8 days, 6 hours. We show in Figure 3 that our 100% networks also achieve similar performance compared to the metrics published by Tor [67]. Our results are plotted with 95% confidence intervals to better understand how well our sampling methods are capable of reproducing the performance characteristics of the true Tor network. We describe how to conduct such a statistical inference in §5 next.

A single network sampled from the true Tor network may not consistently produce perfectly representative results due to the sampling error introduced in the model sampling process (i.e., §3). Similarly, a single simulation may not perfectly represent a sampled network due to the sampling error introduced by the random choices made in the simulator (e.g., guard selection). Multiple samples of each are needed to conduct a statistical inference and understand the error in these sampling processes.

We independently sample $n>0$ Tor networks according to §3.2. The $i$th resulting Tor network is associated with a probability distribution $\widehat{P}_{i}(X)$ which is specific to the $i$th network and the relays that were chosen when generating it. To estimate $\widehat{P}_{i}(X)$, we run $m_{i}>0$ simulations in the $i$th Tor network. During the $j$th simulation in the $i$th network, we sample $\nu_{ij}$ values of $X$ from $\widehat{P}_{i}(X)$ (i.e., we collect $\nu_{ij}$ time to last byte measurements from the simulation). These $\nu_{ij}$ samples form the empirical distribution $\widetilde{E}_{ij}(X)$, and we have $\sum_{i=1}^{n}m_{i}$ such distributions in total (one for each simulation).

Once we have completed the simulations and collected the $\sum_{i=1}^{n}m_{i}$ empirical distributions, we then estimate the inverse distributions $\widehat{F}_{Xi}^{-1}$ and $F_{X}^{-1}$ associated with the sampled network and true probability distributions $\widehat{P}_{i}(X)$ and $P(X)$, respectively.

First, we estimate each $\widehat{F}_{Xi}^{-1}(y)$ at quantile $y$ by taking the mean over the $m_{i}$ empirical distributions from network $i$:

We refer to $\widehat{\mu_{i}}$ as an estimator of $\widehat{F}_{Xi}^{-1}$; when taken over a range of quantiles, it allows us to estimate the cumulative distribution $\widehat{F}_{Xi}(x)=\widehat{P_{i}}(X\leq x)$.

Second, we similarly estimate $F_{X}^{-1}$ over all networks by taking the mean over the $n$ distributions estimated above:

We refer to $\mu$ as an estimator of $F_{X}^{-1}$; when taken over a range of quantiles, it allows us to estimate the cumulative distribution $F_{X}(x)=P(X\leq x)$.

We visualize the process of estimating $F_{X}^{-1}$ in Figure 4 using an example: Figure 4(a) shows $n=3$ synthetic distributions where the upward arrows point to the $\widehat{F}_{Xi}^{-1}$ values from network $i$ at quantile $y=.5$, and Figure 4(b) shows the mean of those values as the estimator $\mu$. The example applies analogously when estimating each $\widehat{F}_{Xi}^{-1}$.

We quantify the precision of our estimator $\mu$ using CIs. To compute the CIs, we first quantify the measurement error associated with the empirical samples. This will often be negligible, but a possible source of nontrivial measurement error is resolution error; that is, if the empirical results are reported to a resolution of $r$ (e.g., 0.01 s), the resolution error for each sample will be $\frac{r}{\sqrt{12}}$, and the resolution error ${\zeta_{i}}$ for the empirical mean $\widehat{\mu_{i}}(y)$ of network $i$ at quantile $y$ is ${\zeta_{i}}=\frac{r}{\sqrt{12m_{i}}}$. Next, we quantify the sampling error associated with the estimates from Equations 1 and 2. The error associated with $\widehat{\mu}_{i}$ for network $i$ at quantile $y$ is:

where $\widehat{\sigma}_{i}(y)=\sqrt{\frac{1}{m_{i}}\sum_{j=1}^{m_{i}}(\widetilde{F}_{Xij}^{-1}(y)-\widehat{\mu_{i}}(y))^{2}+\zeta_{i}^{2}}$ is the standard deviation over the $m_{i}$ empirical values at quantile $y$ (including the measurement error) and $t$ is the $t$-value from the Student’s $t$-distribution at confidence level $\alpha$ with $m_{i}-1$ degrees of freedom [25, §10.5.1]. $\widehat{\epsilon}_{i}(y)$ accounts for the sampling error and estimated true variance of the underlying distribution at $y$. The error associated with $\mu$ at quantile $y$ is:

where $\sigma(y)=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(\widehat{F}_{Xi}^{-1}(y)-\mu(y))^{2}}$ is the standard deviation over the $n$ estimated inverse distribution values at quantile $y$, and $\delta(y)=\frac{1}{n}\sum_{i=1}^{n}\widehat{\epsilon}_{i}(y)$ is the mean error from $\widehat{\mu}_{i}$ over all $n$ sampled networks. $\epsilon(y)$ accounts for the sampling error introduced in the Tor network model generation and in the simulations. We can then define the CI at quantile $y$ as the interval that contains the true value from the inverse distribution $F_{X}^{-1}(y)$ with probability $\alpha$:

The width of the interval is $2\cdot\epsilon(y)$, which we visualize at $y=.5$ with the downward arrows and over all quantiles with the shaded region in Figure 4(b).

Recall that we collect $\nu_{ij}$ empirical samples of the random variable $X$ from simulation $j$ in network $i$. If we increase $\nu_{ij}$ (e.g., by running the simulation for a longer period of time), this will result in a “tighter” empirical distribution $\widetilde{E}_{ij}(X)$ that will more closely resemble the probability distribution $\widehat{P}_{i}(X)$. However, from Equation 1 we can see that $\widetilde{E}_{ij}(X)$ only contributes a single value to the computation of $\widehat{\mu}_{i}$ for each quantile. Therefore, once we have enough samples so that $\widetilde{E}_{ij}(X)$ reasonably approximates $\widehat{P}_{i}(X)$, it is more useful to run new simulations than to gather additional samples from the same simulation.

Additional simulations in network $i$ will provide us with additional empirical distributions $\widetilde{E}_{i*}(X)$, which will enable us to obtain a better estimate of $\widehat{P}_{i}(X)$. Moreover, it will also increase the precision of the CI by reducing $\widehat{\epsilon}_{i}$ in Equation 3: increasing the number of $\widetilde{E}_{i*}(X)$ values at each quantile will decrease the standard deviation $\widehat{\sigma}_{i}$ (if the values are normally distributed) and the $t$-value (by increasing the number of degrees of freedom) while increasing the square root component (in the denominator of $\widehat{\epsilon}_{i}$).

Additional simulations in independently sampled Tor networks will provide us with additional estimated $\widehat{P}_{i}(X)$ distributions, which will enable us to obtain a better estimate of $P(X)$. Similarly as above, additional $\widehat{P}_{i}(X)$ estimates will increase CI precision by reducing $\epsilon$ in Equation 4: the standard deviation $\sigma$ and the $t$-value will decrease while the square root component will increase.

To give a concrete example, suppose $\sigma$ is normally distributed according to $\mathcal{N}(1,1)$. The width of the resulting CI for each number of sampled networks $n\in[2,100]$ at quantiles $y\in\{0.5,0.9,0.99\}$ (i.e., P50, P90, and P99, respectively) is shown in Figure 5. Notice that the y-axis is drawn at log-scale, and shows that the width of the CI can be significantly reduced by more than an order of magnitude after running experiments in even just a small number of sampled networks. Additionally, we can see that the main improvement in confidence results from the first ten or so sampled networks, after which we observe relatively diminishing returns.

Another important factor to consider is the network scale $0<s\leq 1$. Larger scales $s$ (closer to 1) cause the probability distribution $\widehat{P}_{i}(X)$ of each sampled network to cluster more closely around the true probability distribution $P(X)$, while smaller values cause the $\widehat{P}_{i}(X)$ to vary more widely. Larger scales $s$ therefore induce smaller values of $\sigma(y)$ and therefore $\epsilon(y)$. (See §6.3 for a demonstration of this phenomenon.)

While $\epsilon$ includes the error due to sampling a scaled-down Tor network (i.e., §3), the main error that is accounted for in $\widehat{\epsilon}_{i}$ is the sampling error introduced by the choices made in the simulator. If this error is low, running additional simulations in the same network will have a reduced effect. To check the sampling error introduced by Shadow, we ran 9 simulations (3 simulations in each of 3 independently sampled networks of scale $s=0.1$) with unique simulator seeds. Figure 6 shows that the empirical distributions of the 50 KiB download times vary much more widely across sampled Tor networks than they do across simulations in the same network. Although it is ideal to run multiple simulations in each of multiple sampled networks, our results indicate that it may be a better use of resources to run every simulation in an independently sampled network. We believe this to be a reasonable optimization if a lack of available computational resources is a concern.

We have established a methodology for estimating the true distribution of random variables being studied across simulations in multiple Tor networks. Importantly, our methodology includes the computation of CIs that help researchers make statistical arguments about the conclusions they draw from Tor experiments. As we explained above and demonstrated in Figure 5, running simulations in smaller-scale Tor networks or in a smaller number of Tor networks for a particular configuration leads to larger CIs that limit us to drawing weaker conclusions from the results. Unfortunately, previous Tor research that utilizes Tor networks has focused exclusively on single Tor networks while completely ignoring CIs, leading to questionable conclusions (see §2.4).We argue that our methodology is superior to the state-of-the-art methods, and present in §6 a case study demonstrating how to put our methods into practice while conducting Tor research.

We refer to an experiment as a unique pair of network scale $s$ and load $\ell$ configurations, and a simulation as a particular execution of an experiment configuration. We study our hypothesis with a set of 6 experiments; for each experiment, we run multiple simulations in independent Tor networks so that we can quantify the statistical significance of the results following our guidance from §5.

The Tor network scales that a researcher can consider are typically dependent on the amount of RAM to which they have access. Although we were able to run a 100% Tor network for our evaluation in §4, we do not expect that access to a machine with 4 TiB of RAM, as was required to run the simulation, will be common. Because it will be more informative, we focus our study on multiple smaller network scales with more accessible resource requirements while showing the change in confidence that results from running networks of different scales. In particular, our study considers Tor network scales of 1%, 10%, and 30% ($s\in\{0.01,0.1,0.3\}$) of the size of the true Tor network. At each of these network scales, we study the performance effects of 100% and 120% traffic load ($\ell\in\{1.0,1.2\}$) using a process scale factor of $p=0.01$, i.e., each TGen process simulates $1/0.01=100$ Tor users.

Another important consideration in our evaluation is the number $n$ of simulations to run for each experiment. As explained in §5, running too few simulations will result in wider confidence intervals that will limit us to weaker conclusions. The number $n$ of simulations that should be run typically depends on the results and the arguments being made, but in our case we run more than we require to validate our hypothesis in order to demonstrate the effects of varying $n$. As shown in the left part of Table 4, we run a total of 420 simulations across our 6 experiments (three network scales and two load factors) using two machine profiles: one profile included 4$\times$8-core Intel Xeon E5-4627 CPUs running at a max clock speed of 3.3 GHz and 1.25 TiB of RAM; the other included 8$\times$8-core Intel Xeon E5-4650 CPUs running at a max clock speed of 2.7 GHz and 1.5 TiB of RAM.

We run each simulation using an independently sampled Tor network in order to ensure that we produce informative samples following our guidance from §5. Each Tor network is generated following our methodology from §3 using the parameter values described above and Tor network state files from January 2019. The resulting network composition for each scale $s$ is shown in Table 5.

Each simulation was configured to run for 1 simulated hour. The relays bootstrapped a Tor overlay network within the first 5 minutes; all of the TGen clients and servers started their traffic generation process within 10 simulated minutes of the start of each simulation. TGen streams created by Markov clients were set to time out if no bytes were transferred in any contiguous 5 simulated minute period (the default apache client timeout), or if the streams were not complete within an absolute time of 10 simulated minutes. Timeouts for streams created by benchmarking clients were set to 15, 60, and 120 seconds for 50 KiB, 1 MiB, and 5 MiB transfers, respectively.

The time it takes to download a certain number of bytes through Tor (i.e., the time to first/last byte) allows us to assess and compare the overall performance that a Tor client experiences. We measure download times for the performance benchmarking clients throughout the simulations. We present in Figure 7 the time to last byte for 1 MiB file downloads, while noting that we find similar trends for other file download sizes as shown in Figures 9, 10, and 11 in Appendix D. The CDFs are plotted with tail-logarithmic y-axes in order to highlight the long tail of network performance as is typically used as an indication of usability.

Figure 7(a) shows the result of our statistical analysis from §5 when using a network scale of 1% ($s=0.01$). Against our expectation, our estimates of the true CDFs (i.e., the solid lines) indicate that the time to download 1 MiB files actually decreased after we increased the traffic load by 20%. However, notice the extent to which the confidence intervals overlap: for example, the width of the region of overlap of the $\ell=1.0$ and $\ell=1.2$ CIs is about 20 seconds at P90 (i.e., at $x\in[8,28]$ seconds) when $n=10$, and is about 3 seconds at P90 (i.e., at $x\in[16.5,19.5]$ seconds) when $n=100$. Importantly, the estimated true CDF for $\ell=1.0$ falls completely within the CIs for $\ell=1.2$ and the estimated true CDF for $\ell=1.2$ falls completely within the CIs for $\ell=1.0$, even when considering $n=100$ simulations for each experiment. Therefore, it is possible that the $x$ position of the true CDFs could actually be swapped compared to what is shown in Figure 7(a). If we had followed previous work and ignored the CIs, it would have been very difficult to notice this statistical possibility. Based on these results alone, we are unable to draw conclusions about our hypothesis at the desired confidence.

Our experiments with the network scale of 10% offer more reliable results. Figure 7(b) shows the extent to which the CIs become narrower as $n$ increases from 5 to 10 to 100. Although there is some overlap in the $\ell=1.0$ and $\ell=1.2$ CIs at some $y<0.9$ values when $n$ is either 5 or 10, we can confidently confirm our hypothesis when $n=100$ because the estimated true CDFs and their CIs are completely distinguishable. Notice that the CI precision at $n=10$ and $n=100$ has increased compared to those from Figure 7(a), because the larger scale network produces more representative empirical samples.

Finally, the results from our experiments with the network scale of 30% reinforce our previous conclusions about our hypothesis. Figure 7(c) shows that the estimated true CDFs and their CIs are completely distinguishable, allowing us to confirm our hypothesis even when $n=5$. However, we notice an interesting phenomenon with the $\ell=1.2$ CIs: the CI for $n=10$ is unexpectedly wider than the CI for $n=5$. This can be explained by the analysis shown in Figure 5: as $n$ approaches 1, the uncertainty in the width of the CI grows rapidly. In our case, the empirical distributions from the first $n=5$ networks that we generated happened to be more closely clustered by chance, but $n=10$ resulted in a more diverse set of sampled networks that produced more varied empirical distributions. Our conclusions happen to be the same both when $n=5$ and when $n=10$, but this may not always be the case (e.g., when the performance differences between two experiments are less pronounced). We offer general conclusions based on our results later in this section.

The client download error rate (i.e., the fraction of failed over attempted downloads) helps us understand how additional traffic load would impact usability. Larger error rates indicate a more congested network and represent a poorer user experience. We measure the number of attempted and failed downloads throughout the simulations, and compute the download error rate across all downloads (independent of file size) for each performance benchmarking client. We present in Figure 8 the download error rate across all benchmarking clients. (Note that another general network assessment metric, Tor network goodput, is shown in Figure 12 in Appendix D.)

Figure 8(a) shows the result of our statistical analysis from §5 when using a network scale of 1% ($s=0.01$). As with the client download time metric, we see overlap in the $\ell=1.0$ and $\ell=1.2$ CIs when $n=10$. Although it appears that the download error rates decrease when adding 20% load (because the range of the $\ell=1.0$ CI is generally to the right of the range of the $\ell=1.2$ CI), we are unable to draw conclusions at the desired confidence when $n=10$. However, the $\ell=1.0$ and $\ell=1.2$ CIs become significantly narrower (and no longer overlap) with $n=100$ simulations, and it becomes clear that adding 20% load increases the error rate.

Our experiments with the network scale of 10% again offer more reliable results. Figure 8(b) shows significant overlap in the $\ell=1.0$ and $\ell=1.2$ CIs when $n=5$ simulations and a very slight overlap in CIs when $n=10$. However, based on the estimated true CDF and CIs when $n=100$, we can again confidently conclude that increasing the traffic load by 20% increases the download error rate because the CIs are clearly distinguishable. Notice that the CI precision for $\ell=1.2$ compared to the CI precision for $\ell=1.0$ offers an additional insight into the results: the error rate is more highly varied when $\ell=1.2$, indicating that the user experience is much less consistent than it is when $\ell=1.0$.

Finally, the results from our experiments with the network scale of 30% again reinforce our previous conclusions about our hypothesis. Figure 8(c) shows that the estimated true CDFs and their CIs are completely distinguishable, allowing us to confirm our hypothesis even when $n=5$.

We offer some general observations based on the results of our case study. First, our results indicate that it is possible to come to similar conclusions by running experiments in networks of different scales. Generally, fewer simulations will be required to achieve a particular CI precision in networks of larger scale than in networks of smaller scale. The network scale that is appropriate and the precision that is needed will vary and depend heavily on the experiments and metrics being compared and the hypothesis being tested. However, based on our results, we suggest that networks at a scale of at least 10% ($s\geq 0.1$) are used whenever possible, and we strongly recommend that 1% networks be avoided due to the unreliability of the results they generate. Second, some of our results exhibited the phenomenon that increasing the number of simulations $n$ also decreased the CI precision, although the opposite is expected. This behavior is due to random sampling and is more likely to be exhibited for smaller $n$. Along with the analysis from §5, our results lead us to recommend that no fewer than $n=10$ simulations be run for any experiment, independent of the network scale $s$.