Distribution of waiting times in a M/G/1 queue for fixed number of jobs.

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:

$$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.

• The mean is a number. In this context, "distribution of the mean" does not make sense. Distributions only make sense for random variables. – Michael Jun 5 '16 at 9:57
• Assuming you have $N$ jobs there at time 0, and none more arrive, the mean time to empty is $NE[Y]$ (where $E[Y]$ is the average service time). Here I am assuming all service times are identically distributed with mean $E[Y]$. You could potentially define a mean "waiting time" as $\frac{1}{N}\sum_{i=1}^N i E[Y]$, since the first job waits an average of $E[Y]$ to exit, the second waits an average of $2E[Y]$, and so on. – Michael Jun 5 '16 at 9:59
• I added some info. Also I stipulate that not all jobs are arriving at time 0 (which I believe you think is the case according to your second comment). – hbogert Jun 5 '16 at 13:00
• Hmm. What you are asking for is the distribution of the empirical mean seen over the first $N$ arrivals. That sounds way more complicated than verifying a simulation. One way to verify the simulation is to just check that your service times average out to the correct value, inter-arrival times do also, and check the empirical correlation. – Michael Jun 5 '16 at 17:26
• You can also use the iterative formula: If $W_n$ is the delay in the queue of job $n$ (not including service time), then $W_1=0$ (assuming the system is initially empty) and $$W_{n+1} = \max[W_n + \tau_{n} - I_{n+1},0] \quad \forall n \in \{1, 2, 3, \ldots\}$$ where $\tau_n$ is the service time of job $n$ and $I_{n+1}$ is the inter-arrival time between arrivals $n$ and $n+1$. Your simulator should be getting these exact $W_n$ values. – Michael Jun 5 '16 at 17:28

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, ...\}$$
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)$.