When I was working on Cambricon-Q: A Hybrid Architecture for Efficient Training last year, I found that many graduate students do not know how to estimate errors, nor did they realize that error estimation is needed in randomized experiments. The computer architecture community does not pay enough attention to this. When performing experiments, our colleagues often run and measure a program of around ten milliseconds on the real system with the linux time command, and then report the speedup ratio compared with the simulation results. Such measurements are erroneous, and the results often fluctuate within a range; when the speedup is up to hundreds, a small measuring error may result in a completely different value in the final report. In elementary education, the experimental courses will teach us to “take the average of multiple measurements”, but how many times should “multiple times” be? Is 5 times enough? 10 times? 100 times? Or do you keep repeating it until the work is due? We need to understand the scientific method of experimentation: error estimation in repeated random experiments. I’ve found that not all undergraduate probability theory courses teach this, but it’s a must-have skill for a career in scientific research.
I am not a statistician, and I was also awful at courses when I was an undergraduate, so I am definitely not an expert in probability theory. But I used to delve into Monte Carlo simulation as a hobby, and I’m here to show how I did it myself. Statisticians, if there is a mistake, please correct me.
Gaussian Estimation
We have a set of samples , now assume that they are independent and identically distributed (i.i.d.), subject to a Gaussian distribution .
Although there is in the Gaussian distribution parameters, which is what we want, is not directly reflected in the measured samples. To approximate , we can take the arithmetic mean of the measured sample values. We first sum to get , then divide by . We get the arithmetic mean . Note that when a Gaussian random variable is divided by , the variance will be reduced by a factor of accordingly.
Although the above argument first assumes a Gaussian distribution, according to Levy-Lindeberg Central Limit Theorem (CLT), the infinite sum of any random variable tends to follow a Gaussian distribution, as long as those random variables satisfy the following two conditions:
Finite variances;
Independence.
So we can derive the distribution of directly from the CLT, no matter what distribution the samples actually follow, and the distribution is the same Gaussian.
As long as there are enough samples, will converges to .
The Gaussian distribution tells us that we are 99.7% confident that the difference between and the true value we want is no more than three times the standard deviation, i.e., . We say that the confidence interval for the 99.7% confidence level of estimated is .
As the number of experiments increases, the standard deviation of will continue to reduce. Originally we only knew that averaging over multiple measurements would work, and now we also know the rate of convergence: for every 4x more samples, the error will shrink by a factor of 2.
The constants 3 and 99.7% we use are a well known point on the Cumulative Distribution Function (CDF) of the Gaussian distribution, obtained by table look-up. Similar common values are 1.96 and 95%. To calculate the confidence level from the predetermined confidence interval, we use the CDF: , where the CDF only excludes the right-side tail, while the interval we estimate is two-sided, so the confidence level should subtract the left-side tail, i.e., . If you want to calculate the confidence interval from the target confidence level, you need to use the inverse function of CDF (ICDF), also known as PPF.
When programming, the CDF of the Gaussian distribution can be implemented relatively straightforward with the math library function erf. Its inverse function erfinv is not common, as it does not exist in the standard library of C or C++. However, CDF is continuous, monotonic, and derivable. Its derivative function is the Probability Density Function (PDF), so use numerical methods such as Newton’s method, secant, and bisection can all be applied for inversion efficiently.
However, the above estimation is entirely dependent on the CLT, which after all only describes a case in limits. In the actual experiment, we took 5 and 10 experimental results, which is far from the infinite number of experiments assumed. Do we still apply CLT though? This is obviously not scientific.
Student t-distribution
There is a more specialized probability distribution, the Student t-distribution, to define the sum of finitely many samples from random experiments. In addition to the above Gaussian estimation process, the t-distribution also considers the difference between the observed variance and the true variance .
The t-distribution is directly related to the degrees of freedom, which can be understood as in repeated random experiments. When the degrees of freedom are low, the tails of the t-distribution are longer, leading to a more conservative error estimate; as the degrees of freedom approach infinity, the t-distribution converges to a Gaussian distribution. For example, 5 samples, 4 degrees of freedom, and 99.7% confidence level, Gaussian estimated 3 times the standard deviation, while the t-distribution estimated 6.43 times, which is significantly enlarged.
When the number of trials is small, we use t-distribution instead of Gaussian distribution to describe . Especially when relying on error estimates to determine the early termination of repeated experiments, the use of Gaussian estimates is more prone to premature termination (especially when the first two or three samples happen to be similar).
Implementation
In Python, use packages: for Gaussian estimators use scipy.stats.norm.cdf and scipy.stats.norm.ppf, and for Student t-distribution estimators use scipy.stats.t.cdf and scipy.stats.t.ppf.
It is more difficult to implement in C++, the STL library is not fully functional for statistics, and it is inconvenient to introduce external libraries. The code below gives a C++ implementation of the two error estimators. The content between the two comment bars can be extracted as a header file for use, and the following part shows a sample program: Repeat experiments on until there is a 95% confidence that the true value will be included within plus or minus 10% of the estimated value.
/* error estimator for repeated randomized experiments. every scientists should know how to do this... requires c++17 or above. Copyright (c) 2021,2022 zhaoyongwei<zhaoyongwei@ict.ac.cn> IPRC, ICT CAS. Visit https://yongwei.site/error-estimator Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */
// ==================================================================== // contents before the DEMO rule is a header. // ==================================================================== // #pragma once // <-- don't forget opening this in a header. // compilers don't like this in the main file.
#define __STDCPP_WANT_MATH_SPEC_FUNCS__ 1 // for std::beta #include<cmath> #include<iterator> #include<tuple>
namespace errest {
// gaussian distribution. works as an approximation in most cases // due to CLT, but is over-optimistic, especially when samples are // few. template<classRealType =double> structgaussian { static RealType critical_score(RealType cl, size_t df){
// due to CLT, mean of samples approaches a normal distribution. // `erf` is used to estimate confidence level. // now we set a confidence level and try obtain corresponding // confidence interval, thus we need the inverse of `erf`, // i.e. `erfinv`. using std::log; using std::sqrt; using std::erf; using std::exp;
// solve erfinv by newton-raphson since c++stl lacks it. auto erfinv = [](RealType x) { constexpr RealType k = 0.88622692545275801; // 0.5 * sqrt(pi) RealType step, y = [](RealType x) { // a rough first approx. RealType sign = (x < 0) ? -1.0 : 1.0; x = (1 - x)*(1 + x); RealType lnx = log(x); RealType t1 = RealType(4.330746750799873)+lnx/2; RealType t2 = 1 / RealType(0.147) * lnx; return sign * sqrt(-t1 + sqrt(t1 * t1 - t2)); }(x); do { step = k * (erf(y) - x) / exp(-y * y); y -= step; } while (step); // converge until zero step. // this maybe too accurate for our usecase? // add a threshold or cache if too slow. return y; }; // we hope that the "truth" E[mean] lands in the interval // [mean-x*sigma, mean+x*sigma] with a probability of // conf_level = erf( x/sqrt(2) ) // conversely, // x = sqrt(2) * erfinv(conf_level) returnsqrt(RealType(2)) * erfinv(cl); // to get interval, mult this multiplier on the standard error of mean. } };
// student's t-distribution. yields an accurate estimation even for few // samples, but is more computational expensive. template<classRealType =double> structstudent { static RealType critical_score(RealType cl, size_t df){ // similar to gaussian, we need to evaluate ppf for critical score. // ppf is the inverse function of cdf, and pdf is the derivative // function of cdf. we find the inverse of cdf by newton-raphson. using std::beta; using std::sqrt; using std::pow; using std::lgamma; using std::exp; using std::log; using std::nextafter;
auto pdf = [](RealType x, RealType df) { returnpow(df / (df + x * x), (df + 1)/2) / (sqrt(df) * beta(df/2, RealType(0.5))); };
auto cdf = [](RealType x, RealType df) { // ibeta, derived from https://codeplea.com/incomplete-beta-function-c // original code is under zlib License, // Copyright (c) 2016, 2017 Lewis Van Winkle // i believe this snippet should be relicensed due to thorough rework. auto ibeta = [](RealType a, RealType b, RealType x){ const RealType lbeta_ab = lgamma(a)+lgamma(b)-lgamma(a+b); const RealType first = exp(log(x)*a+log(1-x)*b-lbeta_ab) / a; RealType f = 2, c = 2, d = 1, step; int m = 0; do { RealType o = -((a+m)*(a+b+m)*x)/((a+2*m)*(a+2*m+1)); m++; RealType e = (m*(b-m)*x)/((a+2*m-1)*(a+2*m)); RealType dt = 1 + o/d; d = 1 + e/dt; RealType ct = 1 + o/c; c = 1 + e/ct; f *= (step = (ct+e)/(dt+e)); } while(step - 1.); return first * (f - 1); }; RealType root = sqrt(x*x+df); returnibeta(df/2, df/2, (x + root)/(2*root)); };
cl = 1 - (1 - cl) / 2; // solve the root of cdf(x) - cl = 0. // first use newton-raphson, generally starting at 0 since cdf is // monotonically increasing. disallow overshoot. when overshoot // is detected, an inclusive interval is also determined, then // switch to bisection for the accurate root. RealType l = 0, r = 0, step; while((step = (cdf(r, df) - cl) / pdf(r, df)) < 0) { l = r; r -= step; } while (r > nextafter(l, r)) { RealType m = (l + r) / 2; (cdf(m, df) - cl < 0 ? l : r) = m; } returnnextafter(l, r); // let cdf(return val) >= 0 always holds. } };
// std::tie(mean, error) = estimator::test(first, last, cl) // // perform error estimation over the range [`first`, `last`). // the ground truth should be included in the interval // [`mean-error`, `mean+error`] with a probability of `cl`. template<class InputIt> static std::tuple<RealType,RealType> test (InputIt first, InputIt last, RealType cl = 0.95){
// sampling has 1 less degrees of freedom. auto n = distance(first, last); auto df = n - 1;
// pairwise sum to reduce rounding error. using std::distance; using std::next; using std::sqrt; auto sum = [](InputIt first, InputIt last, auto&& trans, auto&& sum)->RealType { if (distance(first, last) == 1) { returntrans(*first); } else { auto mid = first + distance(first, last) / 2; returnsum(first, mid, trans, sum) + sum(mid, last, trans, sum); } };
auto trans_mean = [n](RealType x)->RealType { return x/n; }; auto mean = sum(first, last, trans_mean, sum); auto trans_var = [mean, df](RealType x)->RealType { return (x - mean) * (x - mean) / df; }; auto var = sum(first, last, trans_var, sum);
// assuming samples are i.i.d., we estimate `mean` as // // samples[i] ~ N(mean, var) // sum = Sum(samples) // ~ N(n*mean, n*var) // mean = sum/n // ~ N(mean, var/n) // the standard error of mean, namely `sigma`, // sigma = sqrt(var/n). // // due to Levy-Lindeberg CLT, above estimation to `mean` // approximately holds for any sequences of random variables when // `m_iterations`->inf, as long as: // // * with finite variance // * independently distributed auto sigma = sqrt(var/n);
// finally, confidence interval is (+-) `critical_score*sigma`. // here critical_score is dependent on the assumed distribution. // we provide `gaussian` and `student` here. // // - `gaussian` works well for many (>30) samples due to CLT, // but are over-optimistic when samples are few; // - `student` is more conservative, depends on number of // samples, but is much more computational expensive. return { mean, Distribution<RealType>::critical_score(cl, df) * sigma }; } };
using gaussian_estimator = estimator<gaussian>; using student_estimator = estimator<student>;
};
// ==================================================================== // DEMO: experiment with RNG, try to obtain its expectation. // ==================================================================== #include<random> #include<vector> #include<iostream> #include<iomanip>
intmain(int argc, char** argv){
// the "truth" hides behind a normal distribution N(5,2). // after this demo converged, we expect the mean very close to 5. // we use student's t-distribution here, which works better for // small sets. // if you change the distribution to gaussian, you will likely // see an increased chance to fail the challenge. std::random_device rd; std::mt19937 gen(rd()); std::normal_distribution<double> dist{5,2}; // <-- the "truth"
while (1) { // draw 5 more samples then test. std::cout << "Single Experiment Results: "; for (size_t i = 0; i < iteration; i++) { double sample = dist(gen); // <-- perform a single experiment. samples.push_back(sample); if (i < 10) std::cout << sample << " "; elseif (i < 11) std::cout << "..."; } std::cout << std::endl;
// test the confidence interval with our estimator. std::tie(mean, error) = errest::student_estimator::test(samples.begin(), samples.end(), confidence_level); // or with structured binding declaration [c++17]: // auto [mean, error] = errest::student_estimator::test(samples.begin(), samples.end(), confidence_level);