5
$\begingroup$

I run 100 simulations with CI at 99% and different seeds. For each simulation I created a report that contains different performance indexes (means) with their CI. However, I noticed that one index over 5 simulations didn't contain the expected mean. Maybe I'm wrong but when I choose 99% confidence level does it mean that 1 simulation over 100 doesn't contain expected mean?

Edit

I can't upload all source code because is too long in fact I realized a discrete-event simulator. So I can upload a simple code to show my doubt. In essence for every simulation (100 altogether) this code samples an exponential distribution (or another it is not important) with mean 2 until precision (half-width CI/average) is 3% or below. For this part of algorithm I adapted code from chapter 8 of Leemis-Park. Then it increments an accumulator if CI is out from expected mean.
import umontreal.ssj.probdist.StudentDist;
import umontreal.ssj.randvar.ExponentialGen;
import umontreal.ssj.randvar.RandomVariateGen;
import umontreal.ssj.rng.MRG32k3a;

import java.util.concurrent.ThreadLocalRandom;

public class Estimate {
    public static void main(String[] args) {
        long outlier = 0;
        for (int simulation = 0; simulation < 100; simulation++) {
            int numObs = 0;
            double level = .99;
            double currentAvg = 0, currentSumVariance = 0;
            long[] seed = new long[6];
            for (int i = 0; i < 3; i++) seed[i] = ThreadLocalRandom.current().nextLong(4294967087L);
            for (int i = 3; i < 6; i++) seed[i] = ThreadLocalRandom.current().nextLong(4294944443L);
            MRG32k3a.setPackageSeed(seed);
            RandomVariateGen randVar = new ExponentialGen(new MRG32k3a(), .5);

            do {
                double diff = randVar.nextDouble() - currentAvg;
                currentSumVariance += diff * diff * numObs / (double) ++numObs;
                currentAvg += diff / numObs;
            } while (numObs < 2 || StudentDist.inverseF(numObs - 1, .5 * (level + 1)) * Math.sqrt(currentSumVariance / (numObs - 1)) > (.03 * currentAvg) * Math.sqrt(numObs));
            double t = StudentDist.inverseF(numObs - 1, .5 * (level + 1));
            double ci = t * Math.sqrt((currentSumVariance / (numObs - 1)) / numObs);

            System.out.println(currentAvg + " ± " + ci);
            System.out.println("Precision (accuracy) = " + (ci / currentAvg * 100) + "%");
            System.out.println("numObs = " + numObs);
            System.out.println("expected mean = " + randVar.getDistribution().getMean());
            System.out.println("===========");
            if (!(currentAvg - ci <= randVar.getDistribution().getMean() && randVar.getDistribution().getMean() <= currentAvg + ci))
                outlier++;

        }
        System.out.println("outlier = " + outlier);
    }
}

Another thing I used a library for RNG and random variate SSJ by Pierre L'Ecuyer.

$\endgroup$
5
  • 2
    $\begingroup$ Maybe you made a mistake? Show your code? 100 is a small number of simulations for assessing an event that should happen 1% of the time. I'd do 10k (or more) simulations. $\endgroup$ Commented Jul 1, 2022 at 19:34
  • 1
    $\begingroup$ I'd be curious if it stabilizes upon running more simulations, perhaps 1000 or 10000. Even if not, it might just be that your method for calculating the confidence interval is known to be somewhat anti-conservative. // Be careful calling something an "expected value", as that phrase has a specific, technical definition in statistics that I don't think you mean (if the community will pardon my pun). $\endgroup$
    – Dave
    Commented Jul 1, 2022 at 19:36
  • 2
    $\begingroup$ Could you clarify what the event "one index over 5 simulations didn't contain the expected value" might be?? $\endgroup$
    – whuber
    Commented Jul 1, 2022 at 19:42
  • 1
    $\begingroup$ Is the "expected mean" an estimate or is it the true population mean? If it is the former then the expected coverage does not apply to it. $\endgroup$ Commented Jul 1, 2022 at 22:22
  • $\begingroup$ @MichaelLew "expected mean" is not an estimate but it is true population mean. $\endgroup$
    – miticollo
    Commented Jul 1, 2022 at 22:51

2 Answers 2

8
$\begingroup$

First, I want to reinforce two caveats from comments: (a) Per @Dave: $n = 100$ may be too small to see a clear pattern. If you have one defective part among 99 good parts in a sample of size 100, then the number $X$ of bad parts you get in 100 individual draws with replacement would be $X \sim\mathsf{Binom}(100, .01).$ Thus the number of bad parts in one experiment with 100 draws might be that I never see the bad part.

set.seed(1234)
x = rbinom(1, 100, .01)
x
[1] 0

But in 100,000 such draws of 100 parts, the average number $Y$ of bad parts seen in $100\,000$ draws is $1.$ Then I might see a number near to $1.$

set.seed(1235)
y = mean(rbinom(10^5, 100, .01)); y
[1] 1.00104 

Second, per @MichaelLew: your goal is to get a confidence interval for the population mean $\mu,$ not for its estimate $\bar X.$

Finally, to look at a particularly simple case, suppose you are taking samples of size $n = 100$ from the normal population $\mathsf{Norm}(\mu, \sigma=1).$ Then, the z confidence interval for $\mu$ based on $n = 100, \sigma=10$ is of the form $\bar X \pm 1.96\sigma/\sqrt{n}$ or $\bar X \pm 1.96.$

Then you would hope to find that the CI covers (includes) the population mean $\mu = 85$ in about 95% of the $100\,000$ experiments, as illustrated by the simulation (using R) below.

set.seed(2022)
mu = 85
a = replicate(10^5, mean(rnorm(100, mu, 10)))
summary(a)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  80.29   84.33   85.00   85.00   85.67   89.46 
 mean((a - 1.96 < 85) & (a + 1.96 > 85))
 [1] 0.95092   # aprx. coverage prob of 95% z CI
$\endgroup$
1
  • $\begingroup$ Isn't the OP actually asking about seeing 5 bad parts in 5 draws, to use your analogy? outlier = bad part $\endgroup$
    – dipetkov
    Commented Jul 8, 2022 at 6:38
1
$\begingroup$

Caveats about the correct terminology aside, you have implemented Algorithm 8.2.3, pp. 364 in "Discrete-Event Simulation" incorrectly. The relevant line is this while predicate:

while ((n < 40) or (t * sqrt(v / n) > w * sqrt(n - 1)))

where w is the required margin of error (half-width of the confidence interval).

In your implementation, you compare the current margin of error to (a function of) the current average:

t * sqrt(currentSumVariance / (n - 1)) > c * currentAvg * sqrt(n)

The details about the constants are not important (when do we use n-1 vs n?) What matters is that you have currentSumVariance on the left side and currentAvg on the right side of the inequality. So both sides of the comparison are stochastic. currentAvg is the running sample mean, so effectively, the target precision varies from iteration to iteration, esp. at the start of the simulation. Most likely your simulations happened to exit the while loop too early due to the unintended randomness in the required precision / margin of error.

You should keep w fixed.

Aside: This method for computing the sample variance one step at a time, rather than with a big sum over the entire sample, is known as Welford accumulator.

To verify, I coded the algorithm in python. (For reproducibility, I fixed the seed.) There are 13 outliers in 1,000 simulations; this agrees well with the significance level α = 0.01.

[1] L. M. Leemis and S. K. Park. Discrete-Event Simulation: A First Course (2006)

import numpy as np
import scipy.stats as stats

np.random.seed(seed=20220608)

def rexp(rate):
    return stats.expon.rvs(scale=1 / rate, size=1)[0]


def qt(p, df):
    return stats.t.ppf(1 - alpha / 2, df - 1)


def t_margin_of_error(n, alpha, sigma2):
    return qt(1 - alpha / 2, n - 1) * np.sqrt(sigma2 / (n - 1))


def run_simulation(alpha, delta, rate, verbose = False):

    currentAvg = 0.0  # estimate
    currentSumVariance = 0.0  # residual sum of squares

    n = 0

    while n < 100 or t_margin_of_error(n, alpha, currentSumVariance / n) > delta:

        diff = rexp(rate) - currentAvg

        currentAvg += diff / (n + 1)
        currentSumVariance += diff * diff * n / (n + 1)

        n += 1

    half_width = t_margin_of_error(n, alpha, currentSumVariance / n)

    if verbose:
        print(f"CI = {currentAvg:.3f} ± {half_width:.3f}")
        print(f"Precision = {half_width / currentAvg * 100:.2f}%")
        print(f"Iterations = {n}")
        print(f"Population mean = {1/rate:.3f}")
        print("===========")

    return np.abs(1 / rate - currentAvg) > half_width


alpha = 0.01
delta = 0.3
rate = 0.5

outliers = 0
draws = 1000

for _ in range(draws):
    outliers += run_simulation(alpha, delta, rate)

outliers, draws
#> (13, 1000)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.