Is there a formula to calculate the distribution of the mean waiting time for a server which processes a workload of $N$ jobs, which behaves like a M/G/1 queue?
I would like to be able to answer the following: what is the chance that for a workload of $N$ items, the mean waiting time is above X. If taking the distribution of the mean is too troublesome, can we relax the problem to give the distribution of the waiting times for $N$ jobs (although I guess this does not relax the problem that much).
I hope the question is clear.
Update
Based on the comments I see that some things are unclear. I hope context of my problem will make it more clear.
I have made a simulator for a server which takes jobs as input with inter-arrival times which are exponentially distributed. The job service times are according to a general distribution. Basically I use random samples to mimic a M/G/1. For large amounts of jobs I see the mean waiting time converge to the value which corresponds to the following formula for the expected waiting time[1]:
$$ E[W] = \frac{(\lambda * E[\tau^2])}{ (2*(1- \rho))}$$
Where $E[\tau^2]$ is the second moment of the service time.
However when I have a smaller amount of jobs the mean waiting time can vary a lot between runs. Is there a formula which grasps this variance as function of the service-time, the arrival rate, and most importantly, the amount of jobs. Basically this would allow me to verify if my simulator is correct.
Here are some thoughts on unbiased estimation of delay via simulating over renewal periods.
Unbiased estimator
Suppose we have an M/G/1 queue that is in initially empty. Define: \begin{align} W_i &= \mbox{queueing delay of job $i$}\\ \tau_i &= \mbox{service time of job $i$} \\ Y_i &= W_i + \tau_i = \mbox{total delay of job $i$} \\ N(t) &= \mbox{Number of jobs in system (queue plus server) at time instant $t$} \end{align}
Define renewal times $\{t_0=0, t_1, t_2, ...\}$, where $t_i$ is the time when busy period $i$ ends. Define $T_i = t_i-t_{i-1}$ as the duration of renewal frame $i$. Then $\{T_i\}_{i=1}^{\infty}$ are independent and identically distributed (i.i.d.) and we have $$ E[T_i] = \frac{1}{\lambda} + E[busy] $$ where $E[busy]$ for an M/G/1 queue can be easily calculated. Renewal theory says that, with prob 1:
$$ \overline{N} = \frac{E[\int_0^{T_1} N(t)dt]}{E[T_1]} $$ Also, Little's theorem says $\lambda \overline{Y} = \overline{N}$. Further, the standard proof technique of Little's theorem would also show that $$ \int_0^{T_1} N(t)dt = \sum_{i=1}^{M_1} Y_i $$ where $M_1$ is the random number of jobs in the first renewal frame. Hence, if $M_k$ is the random number of jobs in the first $k$ renewal frames, we have: $$ \sum_{i=1}^{M_k} Y_i = \int_0^{t_k} N(t)dt \implies E\left[\sum_{i=1}^{M_k}Y_i\right] = kE\left[\int_0^{T_1}N(t)dt\right] = k E[T_1]\overline{N}$$ Thus, define $R_k$ via: $$ R_k = \frac{\sum_{i=1}^{M_k} Y_i}{\lambda E[T_1]k} $$ Then $R_k$ is an unbiased estimator of $\overline{Y}$, meaning that for all positive integers $k$ we get $E[R_k] = \overline{Y}$. This is nice because it means that a finite simulation produces a random variable that has the exact average delay value. In particular: $$ E[R_k] = \overline{Y} = E[\tau]+\frac{\lambda E[\tau^2]}{2(1-\rho)} \quad \forall k \in \{1, 2, 3, ...\} $$
Variance calc
Since $\sum_{i=1}^{M_k}Y_i$ is the sum of times over $k$ independent renewal frames, the variance of $R_k$ is: $$ Var(R_k) = \frac{Var\left(\sum_{i=1}^{M_1}Y_i\right)}{k \lambda^2 E[T_1]^2} \quad \forall k \in \{1, 2, 3, ...\} $$ Assuming the numerator is finite, we get $Var(R_k) = O(1/k) \rightarrow 0$ as $k\rightarrow \infty$. With a bit more work, I think one could compute the numerator $Var(\sum_{i=1}^{M_1} Y_i)$.