\documentclass{article}

\setlength{\oddsidemargin}{-0.5in}
\setlength{\evensidemargin}{-0.5in}
\setlength{\topmargin}{0.0in}
\setlength{\headheight}{0in}
\setlength{\headsep}{0in}
\setlength{\textwidth}{7.0in}
\setlength{\textheight}{9.5in}
\setlength{\parindent}{0in}
\setlength{\parskip}{0.05in}
\setlength{\columnseprule}{0.3pt}
\usepackage{fancyvrb}
\usepackage{relsize}
\usepackage{hyperref}

\begin{document}

Name: \_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_

Directions: {\bf Work only on this sheet} (on both sides, if needed); 
do not
turn in any supplementary sheets of paper. There is actually plenty of 
room
for your answers, as long as you organize yourself BEFORE starting 
writing.

{\bf 1.} Consider the variable {\bf \_e.events}.

\begin{itemize}

\item [(a)] (5) Suppose in executing a SimPy program, the first three
invocations of {\bf yield hold} had hold time (i.e. the third operand
for {\bf yield}) of 0.52, 0.18 and 0.98, all executed at time 0.0.
Show the current value of {\bf \_e.events}, making sure to use the
correct Python syntax.

\item [(b)] (10) Suppose SimPy were to be changed to use a calendar queue
instead of is current approach, which essentially involves a linear
linked list.  State the numbers of all lines in the function
{\bf\_e.\_post()} in {\bf Simulation.py} which would be need to be
changed.  The structure {\bf \_e.events} will not change.

\end{itemize}

{\bf 2.} Consider the program {\bf MachRep4.py} in our PLN unit on
advanced SimPy.  

\begin{itemize}

\item [(a)] (10) For each of the class variables {\bf UpRate} and {\bf K}
in {\bf MachineClass}, state the effect of increasing the value of that
variable (while holding the other constant) on the long-run proportion
of up time.  In each case, answer either (i) increase, (ii) decrease or
(iii) could either increase or decrease.

\item [(b)] (10) There is an error in {\bf main()}.  Identify it.

\item [(c)] (5) Suppose we were to have another error in {\bf main()}, in
the form of inclusion of a {\bf yield} statement, say inserted after the
call to {\bf initialize()}.  Give a complete list of the numbers of the
lines which would execute.

\end{itemize}

{\bf 3.} (10) Consider the program {\bf Bus5.py} in our PLN unit on
analysis of simulation output.  Suppose that in addition to estimating
E(W) and P(W $>$ 6.2), we wish to also estimate P(W $>$ 10 $|$ W $>$
6.2).  Show what code to add in order to compute and print out this
estimate and a margin of error for it.  {\bf Be very clear as to where
your new code is to be inserted.}

{\bf 4.} Write two functions for random number generation, both using
the inverse transformation method.  Assume they will be added to the
Python {\bf random} module.

\begin{itemize}

\item [(a)] (10) A function {\bf rand05(self)}, for generating variates from
the density $0.5t^{-0.5}$ on (0,1).

\item [(b)] (10) A function {\bf randpois(self,mu)}, for generating
variates from a Poisson distribution with mean $\mu$.

\end{itemize}

{\bf 5.} (10) Consider the program {\bf QoS.py} in our PLN unit on
advanced SimPy.  Suppose a line

\begin{Verbatim}[fontsize=\relsize{-2}]
G.Mon = Monitor()
\end{Verbatim}

is added to {\bf main()}, with a corresponding line in {\bf G}.  Add
code which uses the monitor to collect data so that we will later be
able to form a confidence interval for the mean number of video packets
in the queue for the channel, via the regenerative method.  {\bf Be very
clear as to where your new code is to be inserted.}

{\bf 6.} (10) Write a function {\bf regen(self)} to be added to the {\bf
Monitor} class, for use in computing regenerative confidence
intervals.  It will return a tuple {\bf (cntr,rad)}, consisting of the
center and radius of the confidence interval.  It is assumed that the
function {\bf observe()} is called at regeneration points.

{\bf 7.} (10) Suppose we are simulating an M/M/1 queue.  This means there
is a single server, and that service and interarrival times have the
distributions U(0,1) and U(c-0.5,c+0.5), respectively.  The question at
hand is whether to scan the event list from the front or the back.  Give
precise conditions on c under which it is better to start from the back.
(If you cannot give precise conditions, then give vague ones.)  Hint:
You may find our formula $Var(R) = E(R^2) - [E(R)[^2$ useful.

{\bf Solutions:}  

{\bf 1.a}

\begin{Verbatim}[fontsize=\relsize{-2}]
{0.52:[iterator address], 0.18:[iterator address], {0.98:[iterator address]}
\end{Verbatim}

{\bf 1.b}  The data structure {\bf \_e.timestamps} will no longer be a
simple linked list.  Thus line 200 must be changed.  No other changes
are necessary.

{\bf 2.a}  If {\bf UpRate} is increased, the machine goes down faster,
hence decreased up time.  If {\bf K} is increased, the repairperson is
summoned earlier, thus increased up time.

{\bf 2.b}  In line 73, 2 should be {\bf MachineClass.R}.

{\bf 2.c}  Having a {\bf yield} in {\bf main()} makes the latter a
generator.  Thus calling it returns an iterator rather than resulting in
the function executing.  The only line to execute would be

\begin{Verbatim}[fontsize=\relsize{-2}]
if __name__ == '__main__':  main()
\end{Verbatim}

{\bf 3.a}  After line 21, insert

\begin{Verbatim}[fontsize=\relsize{-2}]
gt70 = 0
\end{Verbatim}

Replace line 26 by

\begin{Verbatim}[fontsize=\relsize{-2}]
if wait > 6.2:
   gt62 += 1
   if wait > 7: 
      gt70 += 1
\end{Verbatim}

At the end of {\bf main()}, insert

\begin{Verbatim}[fontsize=\relsize{-2}]
p7062 = gt70/float(gt62)
print 'the estimated value of P(W > 7.0 | W > 6.2) is', p7062
me = 1.96*sqrt(p7062*(1-p7062)/gt62)
print 'its margin of error is', me
\end{Verbatim}

Note in particular that one divides by {\bf gt62} instead of by {\bf
nreps}.

{\bf 4.a}  Let H denote the corresponding cdf, which by integration we
see is $t^{0.5}$ on (0,1).  Then its inverse is $G(u) = u^2$ on (0,1).
Thus our function is

\begin{Verbatim}[fontsize=\relsize{-2}]
def rand05(self):
   u = self.uniform(0,1)
   return u*u
\end{Verbatim}

{\bf 4.b}  Use the material in the section titled ``The Inverse
Transformation Method'' in the section titled ``Generating Random
Numbers from Discrete Distributions'' in our PLN on random number
generators.  

\begin{Verbatim}[fontsize=\relsize{-2}]
from math import exp

def randpois(self,mu):
   lamb = 1.0/mu
   u = self.uniform(0,1)
   q = exp(-lamb)
   k = 0
   while True:
      if u <= q: return k
      k += 1
      q *= lamb/k
\end{Verbatim}

{\bf 5.}  We cannot use {\bf PerSmp} here, as it is entirely unsuited
for the regenerative method.

First we must choose a set of regeneration points.  The
simplest set would consist of the times at which the queue length
becomes zero.  Doing the same thing for a queue length of one would not
work, due to the nonexponential nature of the data packet lengths,
though it would work if we imposed the additional condition that there
is currently no data packet in the system.  Then, during each
regeneration cycle we would have code to keep track of the
time-integrated queue length, {\bf TimeInt}.  Up hitting a regeneration
point, we would call

\begin{Verbatim}[fontsize=\relsize{-2}]
G.Mon.observe(TimeInt)
\end{Verbatim}

{\bf 6.}  The monitor will consist of a list of two-element lists
$[a,b]$, where {\bf a} is the simulated time at which {\bf observe()}
is called and {\bf b} is the argument in that call.  Then

\begin{itemize}

\item the i$th$ value of {\bf b} is $C_i$

\item the difference between the i$th$ and (i-1)$st$ values of {\bf a}
is $D_i$

\item $n$ is {\bf len(self)}

\end{itemize}

Thus the code is

\begin{Verbatim}[fontsize=\relsize{-2}]
from math import sqrt

def regen(self):
   csum = 0
   csum2 = 0
   dsum = 0
   dsum2 = 0
   cdsum = 0
   n = len(self)
   for i in range(n):
      si = self[i]
      ci = si[1]  # C_i
      csum += ci
      csum2 += ci**2
      # D_i
      if i == 0:  di = si[0]
      else:  di = si[0] - self[i-1][0] 
      dsum += di
      dsum2 += di**2
      cdsum += ci * di
   cbar = csum/n
   dbar = dsum/n
   gammahat = cbar/dbar
   sx = csum2/n - cbar**2
   sx += gammahat**2*(dsum2/n - dbar**2)
   sx -= 2*gammahat*(cdsum/n - cbar*dbar)
   return (gammahat, 1.96*sqrt(sx/n)/dbar)
\end{Verbatim}

{\bf 7.} Since there will be as many services as arrivals, the
distribution of hold time H is a mixture of two uniforms with range 1,
one with mean 1 and the other with mean c, with 0.5 weighting on each.
So,

For the service times S, we have

\begin{equation}
E(S^2) = Var(S) + [E(S)]^2 = \frac{1}{12} + 0.5^2 = \frac{1}{3}
\end{equation}

Similarly, for the interarrival times I, we have $E(I^2) = \frac{1}{12}
+ c^2$.

\begin{eqnarray}
Var(H) &=& E(H^2) - [E(H)]^2 \\
&=& 0.5 \left (\frac{1}{3} \right) + 0.5 (\frac{1}{12} + c^2) - 
[0.5(0.5)+0.5(c)]^2 \\ 
&=& \frac{1}{6} + \frac{1}{24} + 0.5 c^2 
- 0.25 (0.5+c)^2 \\ 
&=& \frac{10}{48} - \frac{1}{16} - 0.25c + 0.25 c^2 \\
&=& \frac{7}{48} - \frac{1}{4} \cdot c + \frac{1}{4} \cdot c^2
\end{eqnarray}

Thus the cutoff point is found by setting this to 1 and solving for c.
We get

\begin{equation}
c = \frac{12+\sqrt{2112}}{24} \approx 2.42
\end{equation}

In other words, start at the back if c $>$ 2.42.

Some people got partial credit for giving reasoned, intuitive
explanations as to why if c is ``very large,'' we should start at the
back.

\end{document}