# The Expected Total Number of Tardy Jobs¶

We compute the expected total number of tardy jobs by means of the tools in Section 4.3 of [sarin10:_stoch_sched].

import numpy as np
from scipy.stats import norm
from scipy.integrate import quad

class Job(object):
def __init__(self, mean, variance, due):
self.mean = mean
self.variance = variance
self.due = due
self.C = 0 # completion time
self.varC = 0
self.rv = None


The example of Table 4.1 of Sarin. The schedule below eq. 4.28 is contained in s.

job1 = Job(40., 15., 120)
job2 = Job(20., 15., 80)
job3 = Job(60., 20., 90)
job4 = Job(25., 5., 60)

s = [job4, job2, job1, job3]


Now compute the completion times.

j1 = s
j1.C = j1.mean
j1.varC = j1.variance
for j1, j2 in zip(s[:-1], s[1:]):
j2.C = j1.C + j2.mean
j2.varC = j1.varC + j2.variance

# F_{[C[j]]}(d_{[j]})

for j in s:
j.rv = norm(loc = j.C, scale = np.sqrt(j.varC) )


The number of tardy jobs, see eq 4.24 and 4.26:

E = sum( j.rv.sf(j.due) for j in s)
print E

1.00000000165


The variance

V = sum( j.rv.sf(j.due)*j.rv.cdf(j.due) for j in s )


And the final result.

for i in range(len(s)-1):
j1 = s[i]
muX = j1.C
varX = j1.varC
sigmaX = np.sqrt(varX)
for j in range(i+1, len(s)):
j2 = s[j]
muY = sum(k.mean for k in s[i+1:j+1] )
varY = sum(k.variance for k in s[i+1:j+1] )
sigmaY = np.sqrt(varY)
res = quad(lambda x: norm.pdf( (x-muX)/sigmaX )*\
norm.cdf( (x-j2.due+muY)/sigmaY ), j1.due, np.inf)
EUU = res/sigmaX
V += EUU
V -= j1.rv.sf(j1.due)*j2.rv.sf(j2.due)

print V

1.64856166217e-09