Statistics in Engineering
With examples in
MATLAB® and R
Andrew Metcalfe, David Green, Tony Greenfield,
Mahayaudin Mansor, Andrew Smith and Jonathan Tuke.
Chapter 13 solutions to odd numbered exercises
Exercise 13.1
Exercise 13.3
The matrix is in fact tridiagonal and sparse and a truncation to a \(21 \times 21\) matrix is adequate, since the buffer can only
have \(9, 10\) or \(11\) packets after one time unit ... and can have \(0, 1, 2, ... , 20\) packets after ten time units.
There are many ways to do this, and the following MATLAB code shows one way. Note: state \(11\) represents the state with \(10\) packets and state \(1\) represents the state with zero packets.
P = 0.1 * eye(21);
P(1,1) = 0.6;
v = ones(1,20);
P = P + 0.4 * diag(v,1) + 0.5 * diag(v,-1);
P10 = P^10;
MyProb = P10(11,11)
Exercise 13.5
We have that
\( \displaystyle P(z) = \sum_{{n=0}}^{\infty} \pi_{n} z^{n}\)
so that \( \displaystyle \frac{d}{dz} P(z) = \sum_{{n=1}}^{\infty} n \pi_{n} z^{n-1}\)
and by induction we have that \( \displaystyle \frac{d^{k}}{dz^{k}} P(z) = \sum_{{n=k}}^{\infty} n(n-1)...(n-k+1) \pi_{n} z^{n-k}\)
thus yielding \( \displaystyle \frac{d^{k}}{dz^{k}} P(z)_{z=0} = k! \pi_{k} z^{n-k}\)
From above we have
\( \displaystyle \frac{d}{dz} P(z) = \sum_{{n=1}}^{\infty} n \pi_{n} z^{n-1}\)
and hence \( \displaystyle \frac{d}{dz} P(z)_{z=1} = \sum_{{n=1}}^{\infty} n \pi_{n}\),
which by definition of expectation gives \( \displaystyle \frac{d}{dz} P(z)_{z=1} = E[N]\).
From above we have
\( \displaystyle \frac{d^{2}}{dz^{2}} P(z)_{z=1} = \sum_{{n=2}}^{\infty} n(n-1) \pi_{n}\)
and so
\begin{eqnarray}\frac{d^{2}}{dz^{2}} P(z)_{z=1}
& = & \sum_{{n=1}}^{\infty} n^{2} - \sum_{{n=1}}^{\infty} n \pi_{n}\\
& = & E[N^{2}] - E[N] \\
& = & E[N^{2}] - (E[N])^{2} + (E[N])^{2} - E[N]\\
& = & V(N) + (E[N])^{2} - E[N] \end{eqnarray}
Putting it all together we see that
$$
V(N) = \frac{d^{2}}{dz^{2}} P(z)_{z=1} + \frac{d}{dz} P(z)_{z=1} - \left(\frac{d}{dz} P(z)_{z=1}\right)^{2}
$$
Exercise 13.7
If there are two machines where breakdowns occur independently at rate \(1\) per day and there is one
mechanic then we can set up a Q matrix representation of both scenarios. The equilibrium probability
distribution for a finite Q matrix is then just the normalised left eigenvector corresponding to the zero eigenvalue.
This is easily found in either Matlab or R. (Hint: Try the command
[vecs,vals] = eig(Q) in Matlab)
When we have one mechanic and the repair is in
two stages (Erlang repair of rate \(2\) per day): find the fault, and fix the fault, where the independent rate for both stages is \(4\) per day, we have
the following Q matrix representation
$$
\left[ \begin{array}{rrrrr}
-4 & 4 & 0 & 0 & 0 \\
0 & -4 & 4 & 0 & 0 \\
1 & 0 & -5 & 4 & 0 \\
0 & 1 & 0 & -5 & 4 \\
0 & 0 & 2 & 0 &-2
\end{array} \right]
$$
where state
\(0\) corresponds to where both machines are broken and no fault is yet found,
\(1\) corresponds to where both machines are broken and one fault is found,
\(2\) corresponds to where one machine is broken and the fault undiscovered,
\(3\) corresponds to where one machine is broken and the fault is found,
\(4\) corresponds to where both machines are operational.
Here \(\pi_{1} = (0.0610, 0.1098, 0.2439, 0.1951, 0.3902)\)
one stage (exponential repair of rate \(2\) per day) with a fix the fault rate of \(2\) per day, we have the following Q matrix representation
$$
\left[ \begin{array}{rrr}
-4 & 4 & 0 \\
1 & -3 & 2 \\
0 & 2 & -2
\end{array}\right]
$$
where state
\(0\) corresponds to where both machines are broken,
\(1\) corresponds to where one machine is broken,
\(2\) corresponds to where both machines are operational.
Here \(\pi_{2} = (0.1111, 0.4444, 0.4444)\)
To compare the two models, we sum the first two elements and the third and fourth elements of \(\pi_{1}\)
to get \((0.1707, 0.4390, 0.3902)\). Hence in the first (Erlang repair of rate \(2\) per day) scenario, the probability of two machines
being operational is less than in the second (exponential repair of rate \(2\) per day) scenario. The probability of no machines being
operational is also far greater in the first (Erlang repair of rate \(2\) per day) scenario than in the second (exponential repair of rate \(2\) per day) scenario.
The probability of a single machine being operational is more comparable but is marginally less in the first (Erlang repair of rate \(2\) per day) scenario.
Exercise 13.9
Using the code given, change the size to \(1000\).
The first \(100\) data are tainted by a transient effect caused by the system being empty at the start. By ignoring these first \(100\) data we
hope to have achieved a better realisation of covariance stationarity. Another way to achieve this is to run a pilot simulation, calculate the average
loadings and preload the system at the start of your actual simulation runs with these average values.
