Probabilistic Timing Estimates in Scenarios Under Testing Constraints

Probabilistic Timing Estimates in Scenarios Under Testing Constraints Sergi Vilardell†, Francesco Rossi∗, Gabriele Giordana∗, Isabel Serra§, Enrico Mezzetti†, Jaume Abella†, Francisco J. Cazorla† Barcelona Supercomputing Center, Spain† Universitat Autònoma de Barcelona, Spain§ Aiko Space S.r.l, Italy∗ Abstract Measurement-based probabilistic (MBP) methods like Extreme Value Theory (EVT) and the Markov’s Inequality have been exploited to derive probabilistic Worst-Case Execution Time (pWCET) estimates. Usually, the reliability and accuracy of pWCET techniques have been evaluated on medium to large sample sizes, 𝑁=[ 10 3, 10 5] . However, several works increasingly advocate for containing the cost of carrying out the test campaign by reducing the number of executions (i.e. the sample size) required by pWCET analysis. Specific scenarios, for example, impose inherent limitations on the collection of timing measurements due to cost and availability of appropriate testing facilities. In this work, we analyze the impact of small sample sizes on MBP. Our analysis shows that classical EVT models for tail estimation require a threshold that estimates where the tail of the distribution begins. In low sample scenarios, the uncertainty in determining this threshold can compromise the reliability of EVT estimates. We also assess the impact of small samples on RESTK, a time forecast method based on Markov’s Inequality. Our results with synthetic data and representative kernels show that RESTK provides the best trade-off in terms of trustworthiness and tightness for small samples, partly due to not relying on the selection of any threshold, as opposed to EVT. CCS Concepts •Computer systems organization → Embedded systems;• Mathematics of computing → Probability and statistics;• General and reference → Validation;Verification;Reliability. ACM Reference Format: Sergi Vilardell, Francesco Rossi, Gabriele Giordana, Isabel Serra, Enrico Mezzetti, Jaume Abella, and Francisco J. Cazorla. 2025. Probabilistic Timing Estimates in Scenarios Under Testing Constraints. In The 40th ACM/SIGAPP Symposium on Applied Computing (SAC ’25), March 31-April 4, 2025, Catania, Italy. ACM, New York, NY, USA, Article 4, 11 pages. https://doi.org/10.1145/ 3672608.3707895 © 2025 ACM. This is the AAM The definitive Version of Record published in: Sergi Vilardell, Francesco Rossi, Gabriele Giordana, Isabel Serra, Enrico Mezzetti, Jaume Abella, and Francisco J. Cazorla. 2025. Probabilistic Timing Estimates in Scenarios Under Testing Constraints. In Proceedings of the 40th ACM/SIGAPP Symposium on Applied Computing (SAC '25). Association for Computing Machinery, New York, NY, USA, 560–569. https://doi.org/10.1145/3672608.3707895. Plan S Compliance: CC BY 1 Introduction Digitization in different domains is on the rise, with an increasing number of software applications deployed to realize a richer set of functionalities, from entertainment services to medical devices and manufacturing pipelines [ 44 ]. These software applications involve the use of complex AI-based algorithms [ 40 , 41 ] that demand high levels of performance that can only be provided in an efficient manner with complex heterogeneous Common Off The Shelf (COTS) computing platforms based on MPSoCs (Multi Processor System on Chip). A wide set of these applications must abide by firm or soft realtime requirements [ 8 , 18 ], i.e. not necessarily hard. Software timing still represents an important dimension in these systems because, while sporadic budget overruns do not lead to catastrophic consequences, untimely execution may lead to heavily degraded operation or unusable results. As an illustrative example, let us consider a GPS-based tracking function embedded in a smartwatch: an occasional deadline overrun in refreshing the GPS position may affect the precision of the current position but will not cause the system to fail. Clearly, repetitive deadline overruns may eventually lead to a largely inaccurate position and an unusable track recording. Timing estimates are still required to provide enough confidence in the fact that those functionalities will execute most of the times within the allocated budget. Hence, a quality process is required to assess the accuracy of timing estimates during product development. Measurement-based probabilistic (MBP) methods [5, 15, 34, 47] are increasingly adopted to derive trustworthy timing upper-bounds in those scenarios where, the combined complexity of the hardware platform (MPSoC) and software application prevents the effective application of more consolidated deterministic analysis approaches [ 7 , 27 , 40 , 41 , 48 , 52 ]. MBP methods typically deliver estimates in the form of an execution time distribution that, for any given probability 𝑝 , provides an execution time bound (etb). The etb can only be exceeded with a probability 𝑞 such that 𝑞≤𝑝 . Hence, by choosing a sufficiently high etb, we can obtain an arbitrarily low, and hence, affordable, 𝑝 (e.g., comparable to the acceptable failure rates for physical components in the system) even for hard real-time systems, where the risk of overrun must be reduced drastically (e.g., down to [ 10 −12, 10 −15] ). Some works in the literature [6] provide a more complete formal definition of this problem. MBP methods, mainly focused on analyzing high-integrity systems, require the collection of thousands or even hundreds of thousands of execution time measurements samples to derive trustworthy bounds for those extremely low exceedance probabilities [ 5 , 20 , 47 , 53 , 54 ]. However, collecting such large samples may not be always desirable or even practically feasible. The collection of large experimental data can be impeded by the particularly complex functional and operational scenario of the system under analysis; or, more generally, testing implies exceptional effort and costs, so determining when to stop the testing campaign becomes paramount to avoid waste of time and resources [ 16 , 32 , 38 , 39 ]. In this work, we address two fundamental questions: Can the lack of a large sample prevent the application of probabilistic methods altogether? And, can MBP techniques be used to provide bounds for probability ranges of interest for soft and firm real-time applications (e.g., [ 10 −4, 10 −6] )? The rationale behind those questions is that, while extreme target probabilities are the most reasonable target for hard real-time critical systems, they can be overly demanding for firm and soft real-time ones, where sporadic, infrequent deadline misses can be tolerated. For these systems, we might not be interested in the extreme tail of the execution time distribution, but less extreme timing behavior under relatively higher probabilities, as this allows intercepting occasionally (but expected) high execution times. In these cases, we contend it is reasonable to trade some accuracy for efficiency in the determination of probabilistic timing bounds by applying probabilistic methods on lower sample sizes at the cost of producing bounds that are modeling events with higher probabilities. In this paper, we provide a careful analysis of the use of MBP methods to derive execution time bounds for firm and soft realtime applications in the presence of small sample sizes. We build on the consideration that the main focus of timing prediction in those systems is not on determining extremely rare timing events that may not even exist, but rather on intercepting infrequent but expected timing events. From a probabilistic standpoint, this means we are mainly interested in deriving tight timing estimates for events occurring with probabilities in the range [ 10 −4, 10 −6] . In particular, the contributions of this paper are as follows: (1) We develop on the specific problems brought by the use of small samples in the application of EVT. We study the theoretical impact of small samples in EVT, show how to compute the uncertainty associated with the Maximum Likelihood Estimation (MLE) methods typically used in EVT, and propose to capture it using the Delta Method (Section 3). (2) We show that, at theoretical level, the RESTK method based on the Markov’s Inequality, is not fundamentally affected by the use of small samples. Instead, the problems may arise with the application methodology. In that line, we show how derivation of the maximum value of 𝑘 and the derivation of moments requires to be changed for small samples (Section 4). (3) We evaluate both MBP family of methods on a set of extensive synthetic distributions, a representative kernel (matrix multiplication), and two case studies (railway and space). Results on synthetic and real data confirm that methods based on Markov’s Inequality stand out for their trustworthiness and tightness even in low sample scenarios (Section 5). The rest of this paper is organized as follows: Section 2 develops the problem addressed in this work and introduces the most relevant background. Section 3 develops on the main problems that the use of small samples brings to the application of EVT. Likewise, Section 4 presents how small samples affect the application of the Markov’s ... 10-2 10-3 10-5 10-4 10-6 10-7 10-8 10-10 10-9 10-12 10-11 10-13 10-14 10-15 ... Video Playing Mission Critical Low-Frequency Safety Critical High-Frequency Safety Critical Timing Failure Probability Figure 1: Timing failure probability continuum. Reference ranges that in reality depend on the particular domain and application. Inequality. Section 5 quantitatively compares both methods. The paper concludes with Section 6 that summarizes the main related works and Section 7 that presents the main take away messages from this study. 2 Motivation When it comes to risk assessment across different application domains, the majority of safety regulations are based on the As Low As Reasonably Practicable (ALARP) engineering principle, which establishes that risk tolerance shall be assessed against the practicality of the applied risk-reduction techniques. In the context of testing, the benefit of more precise estimations has to be proportioned to the time and effort required to apply the WCET techniques [ 38 ]. For measurement-based techniques this last aspect is ultimately related to the number of tests/runs, a.k.a. sample size. • Costs. Limitations may arise because applications interact in a complex environment with the external world, which may require specific facilities, tools, or devices [ 16 , 39 ], such in, for example, on-the-road testing in automotive. The cost and/or availability of adequate measurement setups or simply the compelling need to reduce the cost and time devoted to timing analysis, in projects that are often conducted under severe time pressure, push towards short test campaigns [ 32 ]. This can be further compounded by difficulties to fully automate the process that makes the collection of measurements a cumbersome process. All these complications have motivated some methods to ascertain whether enough tests have been performed in a timing analysis campaign [38]. • Benefits. Software timing is a non-functional aspect of execution that is not uniformly perceived across industrial domains. For many mainstream applications, timing is mainly recognized as a quality aspect with faster or slower execution being appreciated as good or degraded overall performance of a product. For other applications, the concepts of timing misbehavior and violations are attached to more or less firm requirements on the maximum execution time taken to deliver a result. The impact of a timing budget overrun may vary depending on the application domain and the execution context, ranging from product deterioration and economic loss to environmental disasters and casualties. This elasticity of timing requirements for time-related applications leads to a probability continuum representing the tolerance against system timing failures [31], illustrated in Figure 1. • Entertainment systems (e.g. video and/or audio playing): timing failure rates in the order of [ 10 −2, 10 −4] can be experienced without affecting the user experience. For instance, assuming 25 FPS video playing, a faulty frame every 4 to 400 seconds is often tolerable. • Payloads in mission-critical systems like unmanned exploratory space missions, with instruments collecting and processing data in real-time: the systems operating those instruments can afford some sporadic failures that may degrade the success of the mission, but do not degrade it noticeably if they are sporadic enough. Hence, failure rates in the order of [ 10 −4, 10 −6] are usually affordable (i.e. once every few minutes or few hours for instruments collecting and processing one image at 25 FPS). • Electronic rear-view mirrors in modern cars, replacing real mirrors with cameras and video displays: in this case, video playing is not just for entertainment purposes but for assisting humans in critical operations such as driving. They require even lower timing failure rates as they inherit some (low) integrity requirements. • Lastly, other systems are safety-critical (high integrity requirements), e.g. controlling the operation of an autonomous driving car: the consequences of a misbehavior of the system timing failure rates are usually related to those defined by safety standards for hardware failure rates that, while application-domain specific, usually range in [ 10 −9, 10 −12] per hour. Depending on the frequency of execution of the task under analysis, this might result in lower target probabilities of timing failure per application. For instance, if the application executes 1,000 times an hour (10 3 ) the target failure rate range would be [10−12,10−15]or lower [15]. Overall, different application domains, and applications within the same domain, exhibit different timing requirements. End users must provide quality assurance so that their products behave as expected. The rise of complex MPSoCs has consolidated the use of MBP methods that must provide adequate estimates with the ALARP number of measurements (samples). The latter answers the need to reduce the application cost of MBP methods and, in some cases, to cope with the unavailability of large samples [ 16 , 32 , 38 , 39]. 3 EVT for Small Samples In all application domains, the use of high-quality input data (e.g., high-res images) and AI software requires the adoption of highperformance MPSoCs to provide the necessary levels of performance. As system complexity has increased, and MPSoCs become more popular, MBP methods gained attraction due to their ability to capture and analyze the exceptional level of variability emerging in execution time measurements of AI software running on MPSoCs [7, 40, 41]. The target of the proposed MBP solutions has been deriving etb predicting the probability of exceeding particularly high execution times, often focusing on those values exceeded with probabilities in the range [ 10 −9, 10 −15] . In order to obtain trustworthy predictions for such low exceedance probability ranges, the number of execution time measurements that need to be collected can be in the order of thousands or even hundreds of thousands [5, 20, 47, 54]. We are interested in exploring the application of MBP solutions in a different context considering higher exceedance probability ranges and reduced sample sizes. While there is no theoretical obstacle to the application of MBP to the new setups, altering the application conditions of the methods (i.e. reducing data available for the predictions and targeting higher exceedance probabilities) has clear implications on the final outcome of each MBP method, and needs appropriate analysis and evaluation. Understanding and assessing MBP methods under these new conditions is the main goal of this paper. Classical extreme value statistics [ 51 ] consists of two fundamental approaches: selecting the maximum values for slices of your data, the Block Maxima approach (BM), or selecting a threshold from which the order statistics above it are considered extremes, the Peaks over Threshold (PoT) approach. In this work we focus on the second one, for two reasons. • First, we are interested in modeling the tail of the distribution instead of the maxima. • And second, it has been shown that for small sample sizes, 𝑁< 100, when both methods use the whole sample, the PoT approach makes more efficient use of the data in general [11]. The PoT approach requires addressing the challenge of defining the selection of the threshold that determines when the tail of the distribution starts. The limiting distribution on the tail of a random variable, as described by the fundamental theorem of EVT, is the Generalized Pareto Distribution (GPD) [ 1 ]. Let us define the GPD law as: 𝐺(𝑦;𝜎, 𝜉)=      1−1+𝜉𝑦 𝜎−1 𝜉if 𝜉≠0, 1−exp−𝑦 𝜎if 𝜉=0, (1) where 𝜎 is the scale parameter and 𝜉 is the shape parameter, also known as the extreme value index. This law considers three types of tails depending on the value of the extreme value index. For 𝜉> 0we have heavy tails, for 𝜉= 0we have exponential tails, and for 𝜉< 0we have light tails. The GPD law has the following constraints for light tails, when 𝜉< 0, on the parameters 𝜎> 0, and 𝑦≥ 0when 𝜉≥ 0and 0 ≤𝑦≤𝜎/𝜉 . In this work we focus on the GPD with 𝜉≤ 0, which for 𝜉= 0is an exponential tail, and for 𝜉< 0is the GPD with light tails (GLT). We restrict to 𝜉≤ 0given that the exponential tail is regarded as an upperbound for extreme quantile prediction [28]. In this section, we develop on the impact that small samples have on three different automated tail threshold selection methods and specifically on their performance on extreme value estimation. The selected methods take a different approach for tail threshold estimation, which ensures a wider and more representative scope to our analysis. (1) First, we consider a semi-parametric (SP) model [ 10 ] based on simultaneously fitting a semi-parametric function on the bulk of the distribution, and a GLT or exponential (normally considered in EVT methods for timing estimation) on the tail, successively adjusting the threshold that maximizes the likelihood with the mixture of the two models. (2) Second, we focus on the QQ method [ 4 ], leveraging a metric that minimizes the distance between the fitted GLT or exponential and the estimated quantiles of the bootstrap sample. (3) Finally, we also use the CV plot [ 30 ] method, based on the so-called coefficient of variation and used in previous WCETrelated works [28]. Overall, these three different approaches are expected to exhibit different thresholds and uncertainties on the estimation. Theoretical uncertainty analysis. Producing high-quantile estimations in low sample scenarios can incur a Bias-Variance problem [ 17 ] given that there may be just few points left after selecting a proper threshold. In this work we aim to assess the impact of low samples in EVT models from a theoretical point of view. In the methodologies we show, we can isolate one part of the model uncertainty. The CV, QQ and SP methods, all need to perform first an estimation of the extreme value index 𝜉 , which represents the first source of uncertainty in the models. Therefore, the uncertainty on the estimation of 𝜉 is a lower bound for the uncertainty of the model. For low sample scenarios, if the uncertainty on the estimation of 𝜉 by itself is already high, then the uncertainty will further propagate within the models. 3.1 CV Plot The basis of this method is in the Coefficient of Variation (CV), which is the ratio between the standard deviation and the expected value of the distribution. Taking the sample variance ˆ 𝑉(𝑋) and the sample mean ˆ 𝐸(𝑋)we have the sample CV: bcv =√︃ˆ 𝑉(𝑋) ˆ 𝐸(𝑋).(2) As we are interested in the tail, we can define the residual CV as the CV of the values above a threshold 𝑢 . Let the excess mean be 𝐸(𝑋−𝑢|𝑥>𝑢) and the excess variation be 𝑉(𝑋−𝑢|𝑥>𝑢) then the sample residual CV is: bcv𝑢=√︃ˆ 𝑉(𝑋−𝑢|𝑥>𝑢) ˆ 𝐸(𝑋−𝑢|𝑥>𝑢),(3) The residual CV is almost constant for a sufficiently high threshold [30], and this is the case for the GLT. Interestingly, as we develop below, if bcvgpd and b 𝜉 are maximum likelihood estimators (MLE), and with 𝑉gpd (𝑋−𝑢|𝑥>𝑢)= 𝜎2 (1−ˆ 𝜉)2(1−2ˆ 𝜉) and 𝐸gpd (𝑋−𝑢|𝑥>𝑢)=𝜎 (1−ˆ 𝜉) the residual CV for the GLT can be expressed as: bcvgpd =𝜎 (1−ˆ 𝜉)√︃(1−2ˆ 𝜉)/𝜎 (1−ˆ 𝜉) =1 √︃(1−2ˆ 𝜉) .(4) For high enough thresholds, the residual CV is constant, thus the residual CV is an estimator of the extreme value index 𝜉 . For exponential tails, i.e. 𝜉= 0, the residual CV is 1, which makes it a particularly useful tool for finding exponential thresholds. For 𝜉≠ 0, one should find the threshold from which the residual CV is constant and then compute 𝜉. Application to small samples. For small samples, the stabilization of the residual CV can be challenging. If we consider samples of the order 𝑁= 100, and consider a threshold at the 75% quantile, it leaves us 25 data points, which is enough to estimate the mean and variance of the sample and hence can lead to proper estimates of the residual CV [23]. We explore the asymptotic variance to show a bound on uncertainty related with the sample size. Concretely, we propose to capture its uncertainty using the Delta Method [ 26 ], which can be applied for defining the uncertainty of a statistic that is asymptotically Gaussian. The Delta method is a generalization of the concept of propagation of error, which is widely used in science for accurate representations of complex uncertainties. In this case, we compute the confidence intervals around the residual CV in Equation 4, where the 𝜉 is computed with MLE, therefore making the residual CV an asymptotically normal estimator. Following this argument, we use the general version of the Delta Method: √𝑁[𝑔(ˆ 𝜃) −𝑔(𝜃)] 𝑑 −→ N(0, 𝜎2[𝑔′(𝜃)]2),(5) considering that the derivative of 𝑔 at 𝜃 exists and is not zero. Accordingly, we express the uncertainty of the residual CV as: √𝑁[bcvgpd −cvgpd]𝑑 −→ N(0, 𝜉2[cv′ gpd (𝜉)]2).(6) Finally, we compute the derivative and obtain an expression for the uncertainty around the estimator as a function of the sample size: bcvgpd ≈ N(cvgpd, 𝜉2|(1−2𝜉)−3|/𝑁).(7) Example . In order to isolate the uncertainty around the residual CV as an estimator of 𝜉 , let’s consider a PoT model with a GLT with threshold at quantile 𝑝= 0 . 9and scale parameter 𝜎= 1, and compute the theoretical extreme quantile 𝑝= 1 − 10 −6 . We choose an extreme value index 𝜉1=− 0 . 5and with Equation 7 we can compute the 95% confidence interval for a sample size of 𝑁𝑢= 10 values, which would be the excess sample size with a threshold at 𝑝= 0 . 9with a base sample of 𝑁= 100. Computing all numbers, the confidence interval for the residual CV on 𝜉1=− 0 . 5(with a value of cv1= 0 . 71) is bcv1=[ 0 . 72 , 0 . 69 ] , which is a difference of 4 . 7%. On a similar example, with 𝜉2=− 0 . 25 (with cv2= 0 . 81) has a residual CV of bcv2=[ 0 . 82 , 0 . 80 ] , which also hits a 2% difference. Recall that we fix all parameters and this is the effect on low sample size only on the estimation of the residual CV, which is used to estimate the threshold. If we account for the uncertainty from the estimation of the extreme value index, then the uncertainty around the residual CV dramatically increases. 3.2 QQ plot metric This method proposes a metric based on the QQ plot, which is to be minimized. In particular, the metric is the mean absolute difference between the fitted model and the sample quantiles on the tail. Based on the methodology shown in [ 12 ], we consider both the bootstrapped (QQB) and non-bootstrapped (QQ) versions, the latter been used in [4]. Let 𝑋𝑢=(𝑋−𝑢|𝑥>𝑢) be the excess sample above a threshold 𝑢 , then let 𝑋𝑏 𝑢 be the bth bootstrap sample. In a work of the state of the art [ 12 ], the selected bootstrapped sample size is set to be the same size as the excess sample 𝑋𝑢 given a threshold 𝑢 . Let 𝑄(𝑝, 𝑋𝑢) be the sample quantile function, which in the same work [ 12 ] is based on linear interpolations [ 24 ]. The probabilities for the quantiles checked are based on two parameters which can be chosen 𝑝𝑗 = 𝑗/(𝑚 + 1 ) : 𝑗 = 1 , · · · , 𝑚 . Then, the QQ plot based metric reads as 𝑑𝑏(𝑢)=1 𝑚 𝑚 ∑︁ 𝑗=1 ˆ 𝜎𝑏 𝑢 ˆ 𝜉𝑏 𝑢h(1−𝑝𝑗)ˆ 𝜉𝑏 𝑢−1i−𝑄(𝑝𝑗, 𝑋𝑏 𝑢),(8) where ˆ 𝜎𝑏 𝑢 and ˆ 𝜉𝑏 𝑢 are the bootstrapped scale and shape for the threshold 𝑢 respectively. The non-bootstrapped version is essentially the same methodology but using only the excess sample 𝑋𝑢 once for each threshold. For large samples, this method provides a good metric for threshold selection. It is able to provide a good compromise between Bias and Variance. Note that, for high thresholds, the fitting might not be tight because of sparse data in the tails. Conversely, for low thresholds, the fitting might not be tight given we are not in the tail. Application to small samples. Therefore, the QQ plot based metric is able to compromise and find proper thresholds for either the GLT or exponential tails. However, in a small sample size scenario, the minimum of the metric, which is used to select the threshold value, might be more elusive. Given the small amount of data and high variance between values, the minimum distance may be too close to the tail, providing less certainty of a proper threshold selection. In the next section, we formalize the uncertainty around the quantile estimation produced by the MLE of the extreme value index. 3.3 Semiparametric model The PoT methodology is based on modeling the tail of a distribution with a parametric model, i.e. a GLT, above a certain threshold 𝑢 . The fitting of a parametric model permits the computation of the likelihood of the sample belonging to the selected model. In the case of PoT, the tail model is known because theoretically the GLT is the appropriate model. Otherwise, the model for the bulk can be parametric (e.g. fitting a Weibull distribution), non-parametric (computing the empirical distribution function), or semi-parametric, which is the one we use in this work. A semiparametric model is a model where there can be a specification of an analytical model with fixed parameters, but another component of the model may not be known and be infinite-dimensional. For instance, the Weibull Mixture model is a semiparametric model because while an analytic distribution with finite parameters, the number of components (i.e. the different Weibulls composing the mixture) is neither bounded nor known. The motivation behind combining a tail model with a bulk model, is to select a threshold that maximizes the likelihood of the models on the bulk and the tail at once. The likelihood for a PoT semiparametric model takes the form: 𝐿PoT glt (𝑋)=𝐿bulk (𝑋)(1−𝐻(𝑢|𝛾))𝐿glt (𝑋)(9) =Ö 𝑥𝑖≤𝑢 ℎ(𝑥𝑖|𝛾)Ö 𝑥𝑖>𝑢 1−𝐻(𝑢|𝛾) 𝜎1+𝜉𝑥𝑖−𝑢 𝜎−(1+𝜉) 𝜉, (10) in the case of the GLT, where ℎ(𝑥𝑖|𝛾) is the model for the bulk with a parameter 𝛾to tune. For the exponential case: 𝐿PoT exp (𝑋)=𝐿bulk (𝑋)(1−𝐻(𝑢|𝛾))𝐿exp (𝑋)(11) =Ö 𝑥𝑖≤𝑢 ℎ(𝑥𝑖|𝛾)Ö 𝑥𝑖>𝑢 1−𝐻(𝑢|𝛾) 𝜎exp−𝑥𝑖−𝑢 𝜎.(12) In that way, comparing models with different threshold for the GLT and exponential is possible because the sample size is fixed. Cabras et al. [ 10 ] propose a semiparametric method considering the GLT or exponential distribution for the tail, and using the Lindsay method [ 25 ] to tune the model ℎ(𝑥𝑖|𝛾) for the bulk of the distribution. A semiparametric model can be fitted to the bulk by adjusting a kernel function into multiple windows in the probability distribution function. The number and size of the windows is what makes the model semiparametric, because they are unknown a priori. Lindsay proposes a method to fit the optimal window size and number to adjust the bulk smoothly. Application to small samples. The methodology presented in [ 10 ] produces satisfactory results, however in a small sample scenario it may be difficult to find a robust estimation for both models of the bulk and the tail. This may lead to high uncertainty in the threshold selection. More formally, as part of the analysis in this work, we compute the uncertainty associated with the Maximum Likelihood Estimation (MLE) methods typically used in EVT. The process of fitting tails of a distribution and estimating quantiles outside of the observed sample produces several sources of uncertainty. Uncertainty equally stems from the threshold selection, the determination of scale parameter and extreme value index of the GLT, and the quantile function used to obtain quantile estimations. In this work, we isolate all sources of uncertainty from the estimation of the extreme value index, thus producing a lower bound for the uncertainty associated to using small samples on EVT techniques. We show how even a single source of uncertainty can be significant when dealing with small samples. For the SP and the QQ method, we use the MLE of 𝜉 for GLT. This estimator produces the lowest variance given that it assumes the Cramer-Rao bound [ 42 , 45 ]. In fact, MLE asymptotically converges in law to a normal distribution N(0, 𝜎2)with mean 0and variance 𝜎2as, √𝑁[ˆ 𝜃−𝜃]𝑑 −→ N(0, 𝜎2),(13) where 𝑁 is the sample size, ˆ 𝜃 is the MLE of the true parameter, 𝜃 . Where the Cramer-Rao bound for 𝜃 is: The Cramer-Rao bound reads as: var(𝜃) ≥ 1 𝐼(𝜃)(14) where 𝐼is the Fisher information defined as: 𝐼(𝜃)=−𝐸𝜕2 𝜕𝜃2log(𝑓(𝑋;𝜃))|𝜃(15) In order to derive the Cramer-Rao bound for the particular case of GLT, we apply the MLE of the extreme value index 𝜉 𝜕𝐿glt 𝜕𝜉 =−1 𝜉+1 𝜉2 𝑁 ∑︁ 𝑖=1log1+𝑦𝑖 𝜙=0(16) where 𝑦 and 𝜙 are the sample and scale parameter respectively for the tail. Then, we obtain the estimator, ˆ 𝜉= 𝑁 ∑︁ 𝑖=1 log(1+𝑦𝑖/𝜙).(17) Now, we can compute the Fisher information: 𝜕2𝐿glt 𝜕𝜉2=1 𝜉2−2 𝜉3 𝑁 ∑︁ 𝑖=1 log1+𝑦 𝜙=1 ˆ 𝜉2−2 ˆ 𝜉3 ˆ 𝜉=−1 ˆ 𝜉2(18) Finally, the Cramer-Rao bound is as follows: var(𝜃) ≥ 1 𝐼(𝜃)=1 −(− 1 ˆ 𝜉2) =ˆ 𝜉2(19) This result shows us the lowest variance on the estimator of 𝜉 , which we plug into Equation 13 to obtain its confidence interval: √𝑁[ˆ 𝜉−𝜉]𝑑 −→ N(0, 𝜉2).(20) This expression for the GLT captures how a low sample size can affect the uncertainty around the estimation. Example . Let’s assume a sample size of 𝑁= 100, simulating a threshold at quantile 𝑝= 0 . 9, we have an excess sample size of 𝑁𝑢= 10 to perform MLE on 𝜉 . We compute the 95% confidence interval for 𝑁𝑢= 10, for 𝜉1=− 0 . 5, which gives ˆ 𝜉1∈ [− 0 . 34 ,− 0 . 65 ] ; and 𝜉2=− 0 . 25, which gives ˆ 𝜉2∈ [− 0 . 21 ,− 0 . 29 ] . If we compare these confidence intervals to a bigger sample scenario, with 𝑁= 1000, simulating a threshold at quantile 𝑝= 0 . 9, resulting in an excess sample size of 𝑁𝑢= 100, we observe it yields ˆ 𝜉1∈ [− 0 . 45 ,− 0 . 55 ] and ˆ 𝜉2∈ [− 0 . 23 ,− 0 . 26 ] which are much narrower confidence intervals. To isolate the effect of the uncertainty around 𝜉 on the extreme quantile estimation let us consider a PoT model with a GLT with threshold at quantile 𝑝= 0 . 9and scale parameter 𝜎= 1, and compute the theoretical extreme quantile 𝑝= 1 − 10 −6 . We use the confidence intervals on ˆ 𝜉1 and ˆ 𝜉2 and propagate the error to the quantile estimation 𝑁𝑢= 10. The confidence intervals for quantile 𝑝= 1 − 10 −6 on ˆ 𝜉1 is ˆ 𝑞1∈ [ 5 . 47 , 4 . 39 ] , for a 24% difference; while for ˆ 𝜉2 it is ˆ 𝑞2∈ [ 3 . 87 , 2 . 52 ] , resulting in a 53% difference. With the examples above, we show the minimum variance we would see in these methodologies for low sample scenarios, given that in the examples we are not accounting for the additional variance incurred by the propagation of uncertainty of each methodology. Overall, from the analyses we have presented we conclude that the small sample size has a large impact on the uncertainty of the models by definition. Therefore, there is need to assess the impact of low samples in an empirical way. 4 Markov’s Inequality with Small Sample Sizes 4.1 Introduction to Markov’s Inequality and RESTK Markov’s Inequality (MI) [ 3 ] pertains to the probabilistic concept of concentration inequalities, where the probability of an independent random variable being higher or lower than a certain value 𝑏 can be loweror upper-bounded. In the case of MI, the expression upperbounds the cumulative probability as 𝑃(𝑋≥𝑏) ≤ 𝐸(𝑋) 𝑏,(21) where 𝐸(𝑋) is the expected value and 𝑏 is a value pertaining to the distribution of the random variable 𝑋>0. This particular form of MI is not usable in practice for extreme values due to its pessimism. Vilardell et al. [ 54 ] use the particular case of MI with the power-of-k function (MIK), 𝑃(𝑋≥𝑏) ≤ 𝐸(𝑋𝑘) 𝑏𝑘,(22) where 𝐸(𝑋𝑘) are the theoretical moments of the distribution. Given that, in general, the theoretical moments of a distribution are unknown, the authors adapted the theoretical version of the inequality using the theoretical moments, to the version where the moments are estimated from the sample, 𝑃(𝑋≥𝑏) ≤ 1 𝑁Í𝑁 𝑖𝑥𝑘 𝑖 𝑏𝑘,(23) where Í𝑁 𝑖𝑥𝑘 𝑖/𝑁 is the sample moment estimator [ 43 ]. This adaptation of MIK with the sample moment estimator permits the upperbounding of any positive random variable 𝑋> 0. Unfortunately, tightly upper-bounding extreme values requires increasing the value of 𝑘 , and the sample moment estimator, while being a good estimator for higher order moments, is inconsistent for increasing values of 𝑘 . Vilardell et al. [ 54 ] devise the Restricted k(RESTK) algorithm, where a maximum value of 𝑘 , named max𝑘 is fixed for each target probability to ensure upper-bounding. In RESTK, the central part of the algorithm consists of finding the max𝑘 curve, which indicates the maximum 𝑘 to be used in the sample moment estimator in order not to underestimate. The problem is that for quantiles outside the sample range, it is not possible to assess the underestimation. The key finding is that the relationship between max𝑘 and the logarithm of the probability is linear. Therefore, RESTK can find max𝑘 for the probabilities within the sample size, and then project max𝑘 toward lower probabilities outside the sample with a linear projection. In Algorithm 1, we show the pseudo-code necessary to compute the max𝑘 curve. It builds on a given sample of size 𝑁 and the target probabilities to upper-bound 𝑝target that will be outside of the range of the sample. Then, RESTK selects the probabilities of the quantiles within the sample, 𝑝samp , which should be smaller than 1 /𝑁 for accurate estimates. For each probability within the sample size selected, in Line 3, RESTK uses the quantile estimator suggested in [ 24 ] given that it is approximately unbiased in the median. Then, for each quantile estimated, 𝑞est , RESTK finds the smallest value of 𝑘 that does not underestimate. In order to achieve that, RESTK performs 𝑛bootsims number of bootstrap simulations, with each bootstrap sample, samp𝑏 , of size 𝑛𝑏 . Line 6 shows how RESTK finds the minimum between the quantile estimation made by MIK with the sample moment estimator shown in Equation 23, and the estimated quantile from the sample, |MIK(𝑘, sampb, 𝑝) − 𝑞samp| . The minimum is obtained optimizing 𝑘 for the bootstrap sample 𝑛𝑏 . Finally, the algorithm keeps the smallest 𝑘 from all bootstrap simulations, which will be the max𝑘 for that probability. Once the max𝑘 are Algorithm 1 RESTK algorithm to obtain the max𝑘curve. 1: function RESTK(sample,𝑛b,𝑛bootsims,𝑝samp,𝑝target ) 2: for 𝑝∈𝑝samp do 3: 𝑞samp ←sampleQuantileEstimation(sample, 𝑝) 4: for sim ∈𝑛bootsims do 5: sampb←bootstrap(sample,𝑛b) 6: current𝑘(𝑝) ← findMinOnK|MIK(𝑘, sampb, 𝑝) −𝑞samp | 7: if current𝑘(𝑝)<max𝑘(𝑝)then 8: max𝑘(𝑝)=current𝑘(𝑝) 9: end if 10: end for 11: end for 12: max𝑘(𝑝target ) ← linearAdjust(max𝑘(𝑝samp ), 𝑝target ) 13: return max𝑘(𝑝target ) 14: end function computed for all 𝑝samp , RESTK constructs the max𝑘 curve in Line 12 by computing the linear projection to the target probabilities 𝑝target in the linear form max𝑘(𝑝)=𝛼log(𝑝)+𝛽 , given that the max𝑘 and 𝑝target have been estimated. Then RESTK estimates the max𝑘(𝑝target) using the estimated parameters of the linear adjust, ˆ 𝛼, ˆ 𝛽 and 𝑝target . Finally, RESTK can upper-bound the quantiles with probability 𝑝target by computing MIK with a maximum value of max𝑘(𝑝target). RESTK has been tested for low probabilities between 𝑝=[ 10 −12, 10 −15] , which required sample sizes of 𝑁= 10 6 [ 54 ]. Although the results were satisfactory for extremely high quantiles using big sample sizes, it remains to be understood how to use RESTK in small sample scenarios. 4.2 Impact of Small Samples on RESTK In contrast to EVT, MI uses the whole sample and, therefore, there is no need of filtering out the data using a threshold value. However, the parametrization of MIK with RESTK used on large sample sizes cannot be used in the same way on small ones. max 𝑘 estimation. In RESTK, the key part of the algorithm is finding the max𝑘 curve, determining the maximum 𝑘 to be used in the sample moment estimator in order not to underestimate.The probabilities of the quantiles estimated within the sample were 𝑝samp ={ 10 −3, 10 −4, 10 −5} . Those quantiles were used for the linear projection for target probabilities outside the sample range probabilities (e.g., 𝑝target =[ 10 −12, 10 −15] ). Such projection carries some uncertainty, which exacerbates as we lower the target probability and make it be farther away from the quantiles estimated from the sample. Therefore, the sample size will determine how low the probabilities can be for a robust estimation. In order to find the max𝑘 curve for a small sample size we need to estimate quantiles within the sample range. In the case of 𝑁= 100, for instance, we can choose quantiles at probability 𝑝={ 0 . 95 , 0 . 975 , 0 . 99 } . Additionally, finding the smallest max𝑘 that does not underestimate, in order to prevent pessimism, requires bootstrapping the sample for a robust estimation. We also introduce a change in the bootstrap sample size to 𝑛b= 10 with 𝑛bootsims = 100 bootstrap simulations. However, those small sample sizes can lead to significant variance across simulations. We solve this challenge by making use of the Binomial lower confidence interval (LCI) [ 9 ] to adjust the quantile estimation. The binomial LCI takes into account how many data points are above the quantile selected, and provides a bound for the value of the quantile that corrects the lack of data Figure 2: Comparison for the scaled moment computed from the sample moment estimator for multiple sample sizes for a Weibull distribution of shape 𝛼=8and scale 𝜆=80. points when the sample size is small. Applying the LCI allows us to obtain more consistent quantile estimates from sample to sample, hence avoiding the impact of variance across simulations. Finally, the assessment to do for the adaptation of RESTK is whether or not the sample moment estimator is still accurate for small sample sizes. Moment estimation. When performing estimations on unseen extreme quantiles, the sample size plays a key role in determining how trustworthy the estimation can be. A small sample size is not a huge concern when dealing with summary statistics like the mean, where it can be argued that 𝑁= 30 is enough for a good estimation [ 23 ]. Estimating the moments of a distribution from a sample is a calculation that has similar performance in a small sample scenario. High order moments, i.e. with high value of 𝑘 , are actually hard to estimate even with large amounts of data. MIK needs limits to the value of 𝑘 because the sample moment estimator in Equation 23 is not accurate for high values of 𝑘 . That is why the RESTK methodology is necessary, otherwise, with Equation 23, MIK could safely upper-bound. On the positive side, for low values of 𝑘 , the moments can be estimated satisfactorily with a small sample size. A bigger sample size is useful for extremely high quantile estimation like 𝑝= 10 −15 , but on much lower quantiles, increasing the sample size may not provide with increased tightness. Example . Let’s assume a Weibull distribution for which we compare three sample sizes 𝑁={ 100 , 1000 , 10000 } . We compare the performance of the sample moment estimator for the moments between 𝑘=[ 1 , 30 ] , which at the end-point could be considered high order moments. For each moment 𝑘 , we compute the ratio between the sampled moment, ˆ 𝐸(𝑋𝑘) , and the theoretical moment, that is the scaled moment ˆ 𝐸(𝑋𝑘)/𝐸(𝑋𝑘) , and compute the bias. Figure 2 shows the mean scaled moment, which in fact shows the following: i) the sample moment estimator performs similarly for the two highest sample sizes, 𝑁= 1000 and 𝑁= 10000, hence showing that, for moments of that order, an increased sample size does not provide more information for this estimation. And ii), the case of the smallest sample size 𝑁= 100 shows that, while there is a slight deviation in terms of performance, as the scaled moment for high order moments is around 1 . 02, the performance is satisfactory compared to the higher sample sizes scenarios. Given that the goal of RESTK is to compensate the inconsistency of the sample moment estimator, RESTK can still be used for small sample sizes since the target probability is relatively high (e.g., 10 −6 ), and hence, the value of 𝑘used in Equation 23 needed for such probability is lower. 5 Experimental Evaluation In this Section we compare the adapted version of RESTK introduced in this work and the comprehensive set of EVT instantiations, adopting different threshold selection methods. 5.1 Experimental Setup In the evaluation we include a synthetic representative set of unimodal and mixture parametric distributions (Section 5.2), a representative matrix multiplication kernel (Section 5.3), as well as an industrial case study in the space domain (Section 5.4), all executing on a representative target hardware, namely the NVIDIA AGX Orin [ 35 ]. While multiple experimental setups are considered, they share some common configurations. The sample size is fixed across all experiments in 𝑁 = 100. All samples collected have undergone an independence Ljung-box [ 37 ] test to ensure the applicability for EVT and Markov methodologies. For all experimental scenario we report the bias, or expected value, of each model showing the tightness of the provided bound for the target 10 −6 or 10 −5 exceedance probability. Results are reported for RESTK and EVT under all threshold selection methods addressed in Section 3 (CV, QQ, QQB, SP) with exponential (EXP) and GPD light tail (GLT) distributions. 5.2 Results: Synthetic Distributions Representative Distribution Selection. It is safe to assume that the target program will always terminate in the context of timing analysis of real-time programs, and hence, that a WCET value will always exist upper-bounding all conceivable program runs [ 54 ]. Therefore, the tail of program’s execution time execution is a light distribution that can be upper-bounded by light and exponential tail distributions [ 28 , 36 , 49 ]. Exponential tails can be regarded as pessimistic in some cases [ 54 ], thus providing for a safe upperbound. Heavy tails have the characteristic of decaying slower than exponential tails. Therefore, they are regarded as overly pessimistic for WCET purposes and they are not evaluated in our work. It is also the case that the CDF of the execution time of many programs show several “peaks” usually related to the cache levels, with program’s execution time varying around these peaks. This results in the so-called mixture distributions, which we cover with a solid set of reference multimodal (mixture) distributions in accordance to the state-of-the-art [ 5 , 13 , 22 ] with different tail profiles to increase representativeness: Gaussian, Weibull, and LogNormal variants with weights 𝑤 = { 0 . 8 , 0 . 2 } . The unimodal and mixture Gaussian distributions provide different cases on exponential tails. Similarly the Weibull is used for light tailed scenarios, and finally the Lognormal and Beta distributions are used to be flexible in tail behavior, ranging from light tails to exponential tails depending on the parametrization chosen. We show the analytical formula and parametrization for each distribution, as shown in Table 1. Experimental Results. As the reference ground-truth value used to assess the model tightness we use the theoretically defined quantile functions for each parametric distribution. For the case-studies, the quantile needs to be computed from a bigger reference sample of 𝑁ref = 10 6 runs, which is only used to estimate the ground-truth, to produce an estimation with great confidence. Then we use the base quantile functions of the R programming Table 1: Uniand multi-modal distributions we use. Unimodal distributions Acronym Family Parameters Probability Density Function (G)aussian Gaussian 𝜇=100,𝜎=10 1 𝜎√2𝜋exp −1 2𝑥−𝜇 𝜎2 (W)eibull Weibull 𝛼=4,𝜆=80 𝛼 𝜆𝑥 𝜆𝛼−1exp−𝑥 𝜆𝛼 (B)eta Beta 𝛼=8,𝛽=1/4𝑥𝛼−1(1−𝑥)𝛽−1Γ(𝛼+𝛽) Γ(𝛼)Γ(𝛽) Mixture distributions with weights 𝑤={0.8,0.2} Acronym Family Parameters Probability Density Function (WM)ixtureWeibull 𝜆={100,150},Í2 𝑖=1𝑤𝑖𝛼 𝜆𝑖𝑥 𝜆𝑖𝛼−1exp−𝑥 𝜆𝑖𝛼 𝛼={8,8} (GM)ixture Gaussian 𝜇={100,120},Í2 𝑖=1 𝑤𝑖 𝜎√2𝜋exp−1 2𝑥−𝜇 𝜎2 𝜎={10,10} (LNM)ixtureLognormal 𝜇={5.0,5.1},Í2 𝑖=1𝑤𝑖1 𝑥𝜎𝑖√2𝜋exp−(ln𝑥−𝜇𝑖)2 2𝜎2 𝑖 𝜎={0.1,0.05} language [ 24 ]. We also increase the confidence of the sample quantile estimations by bootstrapping the whole reference sample and compute the expected value of the sample quantile. Unimodal distributions. In Figure 3, we show the bias of each model on the prediction of the quantile with probability 𝑝= 1 − 10 −6 . The G and W distributions can be summarized together as they share the same pattern. The EXP models (CV, QQ, QQB, SP), tend to be safe with the CV and SP being pessimistic and the QQ and QBB being relatively tighter. RESTK is the tightest model on the G case while on the W is similar to the CV. GLT models tend to underestimate in general and changing the GLT model does not incur in significant difference in the bias. The B case shows a different behaviour: the CV is pessimistic due to not being able to find a proper exponential threshold for this sample size, instead a light tail threshold was selected for the CV. All other methods except the SP and RESTK are very close to tightness 1. While these model could be considered good performants in other domains, for real timing analysis they are deemed as optimistic due to not upper-bounding. Mixture distributions. Figure 4 shows the bias for the mixture distributions. The pattern of the biases is similar to the unimodal cases. The EXP models show the same bias pattern, with CV and SP being more pessimistic than QQ and QQB. The GLT show the same behaviour for the GM and WM , but for the LNM the GLT CV and SP are as tight as the QQ and QQB. For the parametric distributions, we have shown that except some particular cases, the GLT will always tend to underestimate in low sample scenarios, and thus are not an appropriate model for WCET analysis. We also showed that RESTK tends to be the tightest model with EXP QQ and EXP QQB being the second tightest models in general. 5.3 Results: Representative Kernel In this section we present experimental results for a representative kernel running on the NVIDIA AGX Orin [ 35 ]. The AGX Orin encompasses and advance MPSoC with complex computing elements, interconnect, and memory controller. All these blocks cause variability in the execution time that align with the need to apply probabilistic methods. We have analyzed edge applications with varying real-time constraints, see Figure 5. Data from different sensors 1 ○ are fed into Figure 3: Fitting for unimodal distributions Figure 4: Fitting for multimodal distributions applications 2 ○ like radar applications, object detection, data fusion, and navigation. Applications build on different types of neural networks (NN) 3 ○ like recurrent and convolutional. As MPSoCs like the Orin allow consolidating several applications, there can be multiple NNs instances running on the CPU and accelerators 4 ○ . CPU is used when accelerators are busy or the working set is reduced and the overheads to offload the kernels (along with its data) to the accelerator is too high [ 2 ]. For instance radar-based object detection that relies on small matrices and LiDAR-based object detection may find accelerators busy running heavier camera-based object detection. Many of these NN heavily build on the ubiquitous matrix multiplication [ 21 ] 5 ○ . It is fair to say that matrix multiplication is the most pervasive operator for many AI applications across different domains. In fact, it has been show that matrix multiplication accounts for a large share of the execution time (between 67% and 98.5%) across deployments [ 14 , 29 ]. For our experiments, we consider two variants of matrix multiplication, a basic one and atiled version which will produce higher variability due to its block-structure functioning. The two left-most charts in Figure 6 show the results for the basic and tiled matrix multiplication, respectively. The bias pattern in the unimodal and mixture cases is also present in the matrix multiplication kernels. EXP CV and SP are the most pessimistic models with QQ and QQB being tighter, although RESTK is the tightest one. All GLT models underestimate. For the tiled matrix multiplication all methodologies except RESTK are either optimistic or too close to tightness 1 to be considered safe for WCET. 5.4 Results: Space Case Study We experimented with a spacecraft autonomous navigation and docking application where an autonomous Guidance, Navigation 1 2 3 4 5 Radar, LiDAR, Camera Prediction, Fusion, Navigation, Radar Convolutional NN, Recurrent NN CPU, GPU, ... Matrix Multiplication Figure 5: From sensors/applications to matrix multiplication. and Control (GNC) system is responsible for acquiring information from the spacecraft asset sensors (cameras, star trackers, inertia measurement units, etc.), assessing the position and altitude of the spacecraft, and issuing a specific maneuver. The autonomous GNC deploys a deep learning model based on convolutional neural networks to perform a docking manoeuvre to an uncooperative target (e.g., space station) on a specific docking site. Among the functionalities covered by the GNC, in our experiments we focused on the pose estimation module. This component takes images from a monocular grayscale camera of the target, and must compute its pose, that is the relative position and rotation (three dimensions each) between the target and the chaser. The pose module results are later used for trajectory planning and actuators. The model data set comprises over 10K satellite distances and configurations, and diverse backgrounds featuring Earth and various lighting conditions. The model is implemented on PyThorch environment and has been trained using a lightweight variant of the Single Shot MultiBox Detector (SSD) [ 55 ], offering a good tradeoff between performance and accuracy. The inference model is executed on the Orin AGX platform. It is worth noting that the application heavily builds on convolutional neural networks whose execution is performed on the GPU. We perform 𝑁ref = 10 6 runs, which took weeks to execute, in order to accurately estimate the reference quantile with probability 𝑝= 10 −5 . For all EVT methodologies and RESTK we used 𝑁= 100 and the same configuration for the Orin. In Figure 6 (right chart), we show the results for the quantile estimations. We observe how the space study shows the same pattern, but in this case EXP QQ and QQB are slightly underestimating. RESTK is the tightest model with EXP SP being a close second. GLT models show no different pattern and incur underestimation also in the space study. 6 Related Works Multiple works have tackled extreme timing estimates from a probabilistic perspective, but the majority of them have not assessed their methodologies for a small sample size. To the best of our knowledge, the works that experimented with a small sample size are the following. The work on MBPTA-CV [ 28 ] suggests a minimum sample size of 𝑁= 100, although their experimentation was done with 𝑁= 500. Silva et al. [ 50 ] use a BM approach with at least 𝑁= 150 samples in the experiments. Finally, Reghenzani et al. [ 46 ] work towards selecting a minimum sample size which has enough statistical power when performing a goodness-of-fit