Using test hypothesis to prove I ran a sufficient number of simulations

Question

I simulated a system using discrete-events simulation, more particularly next-event simulation. The system is M/M/1 (infinite FIFO waiting room). Whenever the system is jobless, the server shutdowns. When a job arrives, the server restarts if it was off, but it cannot serve until the restart is finished. Restart time is defined by an exponential distribution.

• $\lambda$ is the arrival rate (expected inter-arrival is $1/\lambda$)
• $\mu$ is the service rate (expected service time is $1/\mu$)
• $\theta$ is the restart rate (expected restart time is $1/\theta$)

The simulation stops at time $\tau$. Jobs that entered before that time are guaranteed to be served. $\mu > \lambda$ is guaranteed.

I want to observe :

• $\mathbb{E}[B]$, the expected service time
• $\mathbb{E}[S]$, the expected sojourn time (departure - arrival)
• $P_{SETUP}$, the probability that the server is restarting when a client arrives
• $P_{OFF}$, the probability that the server is off when a client arrives

I know the theoretical values: $$\mathbb{E}[S] = \frac{1 / \mu}{1 - \rho} + \frac{1}{\theta}$$ $$P_{OFF} = (1 - \rho) \frac{1/\lambda}{1/\lambda + 1/\theta}$$ $$P_{SETUP} = (1 - \rho) \frac{1/\theta}{1/\lambda + 1/\theta}$$

I've programmed this system. Now, I'm wondering how many simulations is enough so that I can conclude the measures obtained from my simulations are representative of the real system.

I'm wondering if I could draw $n$ samples, then $n_2$ samples, with $n_2 = 2 * n$ and test the null hypothesis $H_0 : v = v_2, H_1 : v \neq v_2$. With $v$ any of the above measures from the first batch of $n$ simulation and $v_2$ the same measure from the second batch of $n_2$ simulations. But I fear that it will only prove that doubling the number of simulations is not significant.

I've also seen that we can do the test with only one sample if we have the population measure (and I do) but this gets me nonsensical results. By example if I do

scipy.stats.ttest_1samp(np.random.exponential(1, 1000000), 1)

I get a statistic and pvalue very different from one run to the other and very rarely less than 0.05 which I think is my purpose.

Answer 1

Ok, I found on how to achieve it and in fact, using

scipy.stats.ttest_1samp(samples, samples_expected_mean)

was indeed the way to go. That's what @ThreeDiag was suggesting. This is a bilateral statistical test using student law.

This gives the pvalue and I only have to assert the pvalue of the test and accept the null hypothesis if $pvalue > \alpha$ or reject it otherwise.

So why didn't I see it sooner as I've already tried that ?

It worked fine and my data was wrong... The problem is that the process I'm measuring wasn't stationary. I had to discard the first half of my raw measures (before aggregating into means) to have something meaningful.