Problem 1

  1. Give the Markov transition matrix for X.

Here we will construct our transition matrix, \(t\).

suppressWarnings(suppressMessages(library(expm)))
suppressWarnings(suppressMessages(library(knitr)))
t <- matrix(c(0.9,0,0,1,0.05,0.85,0,0,0.03,0.09,0.9,0,0.02,0.06,0.1,0),nrow=4)
rownames(t) <- c("from low","from medium","from high","from failed")
colnames(t) <- c("to low","to medium","to high","to failed")
kable(t)
to low to medium to high to failed
from low 0.9 0.05 0.03 0.02
from medium 0.0 0.85 0.09 0.06
from high 0.0 0.00 0.90 0.10
from failed 1.0 0.00 0.00 0.00
  1. A new machine always starts in the low state. What is the probability that the machine is in the failed state three weeks after it is new?

Here we construct our initial state vector, \(i\), and multiply it by the transition matrix to the third power. We will use the matrix exponentiation function in the \(expm\) package. The failure rate will be the last element in the distribution vector. The probability that the machine will be in a failed state three weeks after it is new is 0.0277.

i <- c(1,0,0,0)
i  %*% (t %^% 15)
##         to low to medium   to high  to failed
## [1,] 0.5129639 0.1794735 0.2608973 0.04666541
  1. What is the probability that a machine has at least one failure three weeks after it is new?

With only four possible states and three transitions, we can approach this problem by summing the probabilities of all transitions from \(X_{0} \rightarrow X_{3}\) that include at least one failure.

For the four states, there are 11 transitions from \(X_{0} \rightarrow X_{3}\) that include 1 failure and 1 transition with 2 failures.

\[l \rightarrow f \rightarrow l \rightarrow f = 0.02 \times 1 \times 0.02 = 0.0004\] \[l \rightarrow f \rightarrow l \rightarrow m = 0.02 \times 1 \times 0.05 = 0.001\] \[l \rightarrow f \rightarrow l \rightarrow h = 0.02 \times 1 \times 0.03 = 0.0006\] \[l \rightarrow f \rightarrow l \rightarrow l = 0.02 \times 1 \times 0.9 = 0.018\] \[l \rightarrow h \rightarrow f \rightarrow l = 0.03 \times 0.1 \times 1 = 0.003\] \[l \rightarrow h \rightarrow h \rightarrow f = 0.03 \times 0.9 \times 0.1 = 0.0027\] \[l \rightarrow m \rightarrow f \rightarrow l = 0.05 \times 0.06 \times 1 = 0.003\] \[l \rightarrow m \rightarrow m \rightarrow f = 0.05 \times 0.85 \times 0.06 = 0.0025\] \[l \rightarrow l \rightarrow f \rightarrow l = 0.9 \times 0.02 \times 1 = 0.018\] \[l \rightarrow l \rightarrow l \rightarrow f = 0.9 \times 0.9 \times 0.02 = 0.0162\] \[l \rightarrow l \rightarrow m \rightarrow f = 0.9 \times 0.05 \times 0.06 = 0.0027\] \[l \rightarrow l \rightarrow h \rightarrow f = 0.9 \times 0.03 \times 0.1 = 0.0027\]

The sum of these transitions give us: \(P(N_{Failure} \geq 1) = 0.07085\)

  1. What is the expected number of weeks after a new machine is installed until the first failure occurs?

Here we can take the steady state probabilities, approximated after about 20 transitions, and divide the probability of failure by 1. If we round up the results, we would expect 21 weeks to pass until the first failure occurs.

PIj <- as.vector(i  %*% (t %^% 20))
failure <- PIj[4]
1/failure
## [1] 20.7137
  1. On average, how many weeks per year is the machine working?

Here we will find the probability of the machine working - operating at a low, medium, or high state - for each of the 52 weeks of the year. We will then sum those probabilities to find that it is likely to operate 50 weeks of the year. We get the same result if we use the steady-state probabilities and multiple the probability that the machine is working by the number of weeks.

prob <- 0 
for (week in 1:52){
  PIj <- as.vector(i  %*% (t %^% week))
  prob <- prob + sum(PIj[1:3])
  }
prob
## [1] 49.62941
PIj <- as.vector(i  %*% (t %^% 20))
sum(PIj[1:3])*52
## [1] 49.48958
  1. Each week that the machine is in the low state, a profit of $1000 is realized; each week that the machine is in the medium state, a profit of $500 is realized; each week that the machine is in the high state, a profit of $400 is realized; and the week in which a failure is fixed, a cost of $700 is incurred. What is the long-run average profit per week realized by the machine?

Here we take the steady-state probabilities and multiply it by the vector of the earnings for each state. We will find a long-run average profit of $663 per week.

earnings <- c(1000,500,400,-700)
PIj %*% earnings
##          [,1]
## [1,] 662.6645
  1. A suggestion has been made to change the maintenance policy for the machine. If at the start of a week the machine is in the high state, the machine will be taken out of service and repaired so that at the start of the next week it will again be in the low state. When a repair is made due to the machine being in the high state instead of a failed state, a cost of $600 is incurred. Is this new policy worthwhile?

Here we will modify the transition matrix and earnings vector and see how the long-run profit changes. Average weekly profit the machine is also repaired when it reaches a high state.

t2 <- t # copy transitions matrix
earnings2 <- earnings # copy earnings vector
t2[3,] <- c(1,0,0,0) # modify transitions matrix
earnings2[3] <- -600 # modify earnings vector
PIj2 <- as.vector(i  %*% (t2 %^% 20)) # also reach steady state at 20 transitions
PIj2 %*% earnings2
##          [,1]
## [1,] 770.8894

Problem 2

groups <- c(125,18,20,34)

b <- mean(round(groups/197,3))
round(c(0.5+b/4,(1-b)/4,(1-b)/4,b/4),3)
## [1] 0.563 0.187 0.187 0.063
w <- 0.05
m <- 5000 
burn <- 2000
theta <- numeric(m)

prob <- function(y){
  if (y < 0 || y >= 0.5) return (0)
  return ((2+y)^groups[1])*((1-y)^(groups[2]+groups[3]))*(y^groups[4])
}

u <- runif(m) 
v <- runif(m, -w, w) 
theta[1] <- 0.25
  for (i in 2:m) {
    y <- theta[i-1] + v[i]
    if (u[i] <= min(1,prob(y) / prob(theta[i-1])))
      theta[i] <- y else
        theta[i] <- theta[i-1]
}

index <- 1:m
y1 <- theta[index]
plot(index,y1,type="l",main="",ylab="theta")

xb <- theta[(burn+1):m]
print(mean(xb))
## [1] 0.4800133
mean(xb)
## [1] 0.4800133

For our choice of step-value, we see that the chain is converges to the target distribution after a very short burn-in period and that it mixes will. If we discard the first 200 values, we will get sample mean of the generated chain of approximately 0.48.

Problem 3

Y <- c(4,5,4,1,0,4,3,4,0,6,3,3,4,0,2,6,3,3,5,4,5,3,1,4,4,1,5,5,3,4,2,5,2,2,3,4,2,1,3,2,2,1,1,1,1,3,0,0,1,0,1,1,0,0,3,1,0,3,2,2,0,1,1,1,0,1,0,1,0,0,0,2,1,0,0,0,1,1,0,2,3,3,1,1,2,1,1,1,1,2,4,2,0,0,0,1,4,0,0,0,1,0,0,0,0,0,1,0,0,1,0,1)

fullcondm <- function(m0,lambda0,phi0,yval,lengthY,alpha0,beta0,gamma0,delta0){
  if (m0 > 1){ Lamexp <- sum(yval[1:m0])} else Lamexp <- 0
  if (m0 < lengthY){Phiexp <- sum(yval[(m0+1):lengthY])} else Phiexp <- 0 
  result <- lambda0^(alpha0-1+Lamexp)*exp(-(beta0+m0)*lambda0)*phi0^(gamma0-1+Phiexp)*exp(-(delta0+lengthY-m0)*phi0)
  result
}

n <- length(Y)
iterations <- 5000
mu <- lambda <- phi <- beta <- delta <- numeric(iterations)
mprob <- rep(1/112,112)
mdist <- function(n) {sample(1:112,n,replace=TRUE,prob=mprob)}
mu[1] <- mdist(1)
lambda[1] <- 1
phi[1] <- 1
beta[1] <- 1
delta[1] <- 1
alpha <- 8
g <- 2

for (i in 2:iterations){
  m <- mu[i-1]
  sum1 = sum(Y[1:m])
  if (m+1 > n){sum2 <- sum(Y)} 
    else {sum2 <- sum(Y[(m+1):n])}
  phi[i] <- rgamma(1,sum2+g,delta[i-1]*(n+m)) 
  lambda[i] <- rgamma(1,sum1+alpha,beta[i-1]*(n+m)) 

  for (j in 1:n){
    mprob[j] <- fullcondm(j, lambda[i-1], phi[i-1], Y, n, alpha, beta[i-1], g, delta[i-1])
                }
  beta[i] <- 1 / rgamma(1,alpha,lambda[i-1]+1)
  delta[i] <- 1 / rgamma(1,g,phi[i-1]+1)
  while (beta[i] > 99){
    beta[i] <- 1/rgamma(1,alpha,lambda[i-1]+1)
    }
  while (delta[i] > 99){
    delta[i] <- 1/rgamma(1,g,phi[i-1]+1)
    }
  mprob <- mprob/sum(mprob)
  mu[i] <- mdist(1)
}
mu
##    [1]  83  93 107  97  46  97  97 101  92  97  98  98  97  98  97  98  97
##   [18]  98  97  92  97  97  98  97  97 112  97  49 101  97  97 112  97  41
##   [35]  98  46  97  97  98  97  97  97  97  97  97  98  97  97  98  97  97
##   [52]  97  98  97  97  97  97  97  99  97  46  97  63  97  97  98  97  97
##   [69] 101 110  97  97  98 101  97  97  97  97  97  97  97  97  97  98  46
##   [86]  97  97  97  97  46  97  46  97  46  97  97 112 101 112  97  40  97
##  [103]  44  97  97  97  97  97  98  97  97  97  97  98  97  92  97  98  97
##  [120]  97  97  97  97  97  97  97  97  97  97 112  97  50  97  97  97  98
##  [137]  97  97  97  98  98  97  97 112  97  36  97  97  97  92  98  97 101
##  [154]  97  97  98  97  97  97  97 112  97 111  97  98  97  97  99  98  97
##  [171]  97  97  97  97 103  97  97 111  97  97  46  97   1  97   1  97   5
##  [188]  97   1  99  97  97  97  97  98  98  46  97  92  97  97  97  97  97
##  [205]  97 101  97  97  97  97  98  46  97  91  97  41  97   4  97  97  97
##  [222]  41  97  97  97  97  97  98  97  97  97  97  97  97  97  97  98  97
##  [239]  97  97 102 112  97  40 101  41  97  46  98  42  97  46 101  51  97
##  [256]  46 101  60  92   2  98   1  97   1  46 112  98  41  97  38  97  97
##  [273]  98  97  97  97  97  97  97  97  97  46  97  49  97  46  97  97  98
##  [290]  97  97  98  46  97  97  97  97  97  97 101 112  97  36  97  97  97
##  [307]  97 101  97  97  46 102  96  46  97  39  98  46  97  97  98  46  97
##  [324] 102  98  98  97  97  97  97  97  97  97  97  97  97  97 102  98  97
##  [341]  97  97  97  97  97  97  97  97  97  98  97  92  97 101  97  46  97
##  [358] 104  97  97  97  97 101  97  97  97  97  97  97 101  98  92  97  97
##  [375] 102  98  97  97  97  97  97  99  97  97  97  97 112  97  46  97  97
##  [392]  97  97  98 102  97  97  97  97  97  97 102  97  97  97  97  97  97
##  [409]  97  97  97  97  97  97  97  98  98  97  97  97  97  97  97  97  98
##  [426]  97  97  97  97  97 101  97  46  97   1  97   1 110  64  97  97  97
##  [443]  97  97  97  97 101  97  97  97  97  97  97  98  98  97  97  41 101
##  [460]  97  97  97 112  97  37  97  98  97  46  98  53  97  37  99  93  97
##  [477]  46  97   1 101   1  41   1   1   1   1  97   1  97   1  97 103 102
##  [494]  97  97  97  97  97  97  97 101  98  97  97  97 101  97  97  97  97
##  [511]  97  97  97  97 101  97  97 112  97  40 112  97  36  92  89 102  97
##  [528]  92  97 102  97  97  97  98  97  97  97  99  97  97  98  98  46  97
##  [545]  99  97  97  97  97  97  97  98  97  97  97  97 101  98  97  97 100
##  [562]  97 101  97  97  97 101  97  92  97  92  97  92  97 102  97  97  97
##  [579]  97  97  98  97  98  98  98  97  97  97 101 112  97  41  97  56  97
##  [596]  40 101  97  97  97  99  97  97 102 101  97  97  97  98  98  97  97
##  [613] 112 101 110  97 112 102  56  97  98  97 112 102  38  97  97  98  97
##  [630]  46  97  99 102  92  97  97  97  99 101  97 101  97  97  98  97  97
##  [647]  97  97  97  37  92  97  97  41  97   1  97   1  97   1  98 112  92
##  [664]  36  97   1  97   1  97   1  97   1  99   1  97 107  97  41  97  92
##  [681]  97  93  97 101  97  97  97  98  97  97  97 101  97  97  97  97  97
##  [698]  97  97  98  97 102  97  97 103 101  97  97  97  97  98  97  97  92
##  [715]  97  92  97  97  97  97  99  46  98 100  97  97  97  97  97 101  97
##  [732]  97  98  97  92  97 101 103  97  97  97  97  97  97  97  97  97  98
##  [749]  97  97  97  97  97  98 112  97   3  97  72  97  46  97  40  97  97
##  [766]  97  97  97  97  97  92  98 101  97 112  97  97 112 112  36   7 101
##  [783]  45  97  37  97  41  98  97  97  98  97  97  97  97 101  97  97  98
##  [800]  97  46  97  99  97  97  97  97  97  97  97  97  97  97  97  97  97
##  [817]  97  97  97  97  97  97  97  97  92  97  97  97  97  97  97  97  41
##  [834]  92 101 101  97  97  97  98  46  97 112  97  55  97  97 112  97   3
##  [851]  98   1  97   1  97 101  97  40 104  97  97  92  97  97  97  97  97
##  [868]  98  98  97  98  97  98 112  98  36  97   1  97   1  92   1  97   1
##  [885]  97   1  97  97  97  97  97 101  98  97  97  97  97  41  97   1  97
##  [902]   1  98   1  97   1  97   1  98   1  97  78  97  46  46 102   1  97
##  [919]  92 101  46  97  97  97  97  97  97  97 104  97  92  97  97  97  92
##  [936]  97  98  97  97 101  97  97  97 101  97  97 109  97  97  97  97  97
##  [953]  97  46  97  97 112  97   3  99   1  97   1 112   1   3   1   1 111
##  [970]   1  97   1  97  97  46  97  97  97  97  97  97  98  97  97  97  98
##  [987]  98  97  97  97  97  97  98  97  97  97  97  97 101 101  97  92  97
## [1004]  97  40  92   1  97   1  97   1  97   1  97 112  97  92  97  97 100
## [1021]  46  98 106  97 112  97  41  97 100  97  46  97  63 112  46  45  46
## [1038]  98  98  97  97  97  97  97 101  97  97  92  97  97  97  46  99  97
## [1055]  97  97  98  97  97  97  97  97  97  97  97  97  97  97 102  97  97
## [1072]  97 112  97  18  97   1  97 102 102  97  97 112 112  10  34  98   1
## [1089]  37 110   1  97   1 112   1  28   1  92   1  97   1  97  74  97  46
## [1106]  97  97  97  97  97  98  97  98  97 101 101  92  97  97  98  97  97
## [1123]  97  97  97  99  97  97  97  97  92  97  98  97  97  92  97 112  97
## [1140]  36  97  97  97  46  97  42  99  46  97 112  97  40  97  93  97  46
## [1157]  97 108  97  97  97  97  97  97  98  97  98  97 101  97  97 101  98
## [1174]  97  97  97  98  97  97  98  99  97  46  97  46  46  42  40  97 100
## [1191]  97  97  98  97  41  97 101  98  97  97  97  97  97  97  97  97  97
## [1208]  98  97  97  97  98 110  97  97  97  97  98  97  97  98  97  97  97
## [1225]  97  97  97  97  97  92  97 112  97  61  97  98  97  97  98  98  97
## [1242]  46  97  94  97  46  97   1  97   1  97  93  97  46  97  47  99  46
## [1259]  97  39  98  46  97  47  97  97  97  97  98  97  97  97  97  98  92
## [1276] 101  97  98  46  97  92  97  92  97  97  97 101  97  97  97  97  97
## [1293]  97  97  97  92  97  97  97  97  98 102  97  97  97  97  46  97  97
## [1310]  97  97  97  97  97 112  97  39  97  97  92  99  97  97  98  97  97
## [1327]  97  97  97  97  98  97  97  97  98  97  97  97 102  97  98  97 102
## [1344] 101  97  97  92  97  97 102  98  60  92  97 101 101  46  97   1  97
## [1361]   1  97 109  97  98  97  97  97  97  98  97  97  97  97  98  97  97
## [1378]  97  97  97  97  98 101  92  97  97  92  97  97  97  97  97  98  97
## [1395]  97  97  97  97 112  97  31 101  98 101  97  97  97  97  99 101  97
## [1412]  46  97   5  97   1  97   1  98  92  97  98  97  97  97 101  97  97
## [1429]  97  97  97  97  97  98  97  97  92  97  97  97  92 112  98  92  97
## [1446]  97  97  97  46  97  43  97  45  97  43  97  43  97  99  97  97 100
## [1463]  98  92  46 100 102 101  97  97  97 112  97  37  97  97  97  98  97
## [1480]  97  97  97  97  98  97  92  97  97 101  97  97 110  97  97  97  98
## [1497]  97  97  97 112  98  40  46  97  97  97  92  97  97  98 101  46  97
## [1514]  94  97  97  92  99 107  97  41  97   1 101   1  46   1 107  98  97
## [1531]  46  97  98  97  97  97  97  46  97  99  97  97  97  97  46  97 101
## [1548]  92  97  97  97  97  97  97  97  97  97  97  97  97  97  98  98  97
## [1565]  97  97  97  97 110  97  97  98  97  97  97  97  97  97  97  98  97
## [1582]  97  97  97  97  97  97  97  97  97  97  97  97  97  97  97  97  97
## [1599]  97  97 107 107  97  97  97  97  97  97  98  98  97  98 112  97   3
## [1616]  97   1  46  98  91  97  98 101  97  46 101 100  97  97  97  97  97
## [1633]  97  97  97 101  97  97  97  97  97  97  98  97  97  97 102 101  97
## [1650]  92  97 100  97  97 101  97  97  97  98  97  97  97  97  97  97  97
## [1667]  97 110  97  97  97  97  97  97 112  97  32  97 100  97  97  97  97
## [1684]  97  97  97  97  97  97  97  97  97 101  97  97  97 100  97  92  97
## [1701]  97  97  97  97  97  97  97  97  97 112  99  42  97  47  98  47  97
## [1718]  99  97  97  98  97  97  97  97  97  97  97  97  98  98  97  97  97
## [1735]  97  46  97  61  97  97  97  46  98 101  98  98  97  97  97  97  97
## [1752]  98  97  97 102  97  97  97 112  97  41 101  97  97  97  97  97  97
## [1769]  97  97  97 101  97  97  97  97  97  97  97  97  97  97  97  97  98
## [1786]  97  46  97  38 102  98  46  97  98  97  97  97  97  97  97  97 105
## [1803]  97  97  97 101  97  98  97  97  97  97  97  97  97 101  97  97  97
## [1820] 101  97  97  97  97  97  97  97  97  97  97  97  97  97  98  97  40
## [1837]  92   1  97 101  97  97  97 110  97  92  97  98  97  46  97  97 101
## [1854]  97  46 100  93  97  97  97 112  46  33 101   9  97  84 112  97  41
## [1871]  97  51  97  41  98  99  97  97  92  97  97  97  97  97  97  97  97
## [1888]  97  97  98  97  97  97  97 101  46  92 102  97  46  97  47  97  92
## [1905]  97  98  97  97  97  97  97  97  98  97  97  97  97  98  98  98  98
## [1922]  97  97  97  97  92  97  97  97  97  97  97  46  46   1   1   1   1
## [1939]   1   1   1   4   2   1   1   1   1   1   1   1   1  92  93  92  46
## [1956]  97  93  97  41  97  97  97  97  97  97  97  98  97  98 111  98  97
## [1973]  97  97  97  97  97  97  97  97  97  98  97  97  97  97  97  97 102
## [1990]  98  97  92  97  97  97  97  98  97  97  97 102 110  92  97  97 107
## [2007]  97  97 101  98  46  97  98  98  97  97 102  97  97 102 110  97  97
## [2024]  98 112  97 111  98 112  99  40  97  91  97  43  97   1  92   1  97
## [2041]   1  97 112  97  97  97  97  46  97  42  97  99  97  97  97  97  92
## [2058]  97  99  97  46  97  97  97  97  98  97  97  92 101  97  97  97  98
## [2075]  97  98  97  98  97  97  97  97  97  97  98  98 102  98  97  97  97
## [2092]  97  98  97  97  97  97  97  97  97  97  97  97  97  97  97  98  97
## [2109] 101  97  98  97  99  98  97  97  97  97  99  92  97  46  97  45  97
## [2126]  41  97 102  98  97  98 102  98  97  46  97  92  97  97  97 112  98
## [2143]  34 112  47  40  50  93  46  97  97 102  46  97  98 112  97  46 102
## [2160]  97  97  97 112  97  32  97   1  97   1  99   1 101   1  97   1  98
## [2177] 112 102  32  97   1  97   1  97   1  97   1  97   1  97  92  97  97
## [2194]  46  97  99  97  97  97  97 103  97  92  97  97  97  46  97  98  97
## [2211]  97 112  98  41  97  98  46  97  93  97  98  98  97  97  97  92  98
## [2228] 112  97   2  97   1  97 112  97  46 112  61  37  39  46  97  62  97
## [2245]  99  97  97  97  97  97  97  97  97  97 100  46  97   3  97   1  97
## [2262]   1  98 111  97  97  97  97  98  97  97  97  97  97  97  98  92 110
## [2279] 102  92  46  97  98  97  97  97 110  98 112  97  41 103  42  41  47
## [2296]   1  97   1  97   1 112   1   3   1   1   1   7   1   1 112 107  42
## [2313]  97  15 106 112  97  97  97  97  97  97  97  97  97  97  92 112  97
## [2330]  47  97  97  97  97  97  97  97  98 101  97  46  97   1  97   1  97
## [2347]   1  46   1  97   1 110   1  97   1  97   1  97  97  97  97  97  99
## [2364]  97  97  97 112 112  37  44  92  98  46  97 102  97  46  97  97  97
## [2381]  97  97 100  97  46  98  97  97  97  97 102  97  60  97  91  97  97
## [2398]  97  97  97  97  97  97  97  98  97 102  97  97  97  97  97  97  97
## [2415]  97  97  97  46 101  97  97  97  97  97  92  46 101  42  97  92  97
## [2432]  97  97 101  97  46  39   1  97   1  97  99  97  97  92  97  97  97
## [2449]  98  97  97  97  97  97  97  97  97  98  97  97  97  97  97  97  97
## [2466]  97  98  97  97  97  97  97  98  97  97  97  97  97  97  97  97  97
## [2483]  97 112  97   3  97   1  98   1  97   1  98   1  97 107  97  97  97
## [2500]  97  97  97  97  97  97  97  98 101  97  97  97  98  97 112  46  36
## [2517]  97  97  97  97 101 112  97  41  97  97  97  92  97  97  97  97  97
## [2534]  97  97  92  97  97  97  97  97  97  97  97 101  97  97 101  97  97
## [2551] 112  97   3  97   1  97   1  97  96  97  41  97  97  97  97  97  97
## [2568]  97  97  98  97 112  97 112  97  97  97  97  46  97  39 112  92  46
## [2585]  97  92  47  97  63  97  97 112  97  43  97  92 101  93  97  97  97
## [2602]  97  97  97  97  98  97  97  97  97  97  97  97  97  97  97  97  97
## [2619]  97  97  97  97  97  97  97  97  98  97  99  97  97  97  92  97  97
## [2636] 101  97  97  97 112 112  37  31  98  98  46  98  45 111  46  97  98
## [2653]  97  97  97  97  97  97  97 102  97  97  97  97  97  46  97  97  97
## [2670]  46  97  37  97  97  97  97  97  97 101  98  97 100  97 102  97  97
## [2687]  97  97  97  97  97  98  97  97 101  97  97 101  98  97  97  97  97
## [2704]  97 101  97 112  97  40  41  97  97  98  97  97  97  98  97  97  97
## [2721]  97  97  97  97  97  97  98  97  97  97  46 102  93  92  46  98  98
## [2738]  97  97 101  97  97  97 105  97  97  97  97  97  97  97  97  97  97
## [2755]  97  97  97  97  97  97 112  92  40  97  97  97  46  97  62  97  46
## [2772]  97  92  97  97  46  97  87  97  97  92  97  97  97  97  43  97  97
## [2789]  97  97  97  98 101  97  97  97  97  97 112  97  40  97  97  97  97
## [2806] 101  97  97  97  97  98  97  97  98  97  97  97  97  97  98  97  97
## [2823] 112  97  40  97  49  97  41  97  64  99  97  97  97  97  97  97  98
## [2840]  97  97  97 103  97  97  46  97  42  97  46  97  42  97  97 112  46
## [2857]  28 103  98  97  97 111  97  97  97  97  97  97  97  97  92  97  97
## [2874]  98 112  97  32 101  91  97  97  97  97  97  97  97  97  97  97  97
## [2891]  97  97  97  97  97  46  97  80  97  40  97  93  97  97  97  97  97
## [2908]  99  92  45  98  41  97  46 112  97  42  92  97  97  97  97  97  97
## [2925]  97  97  97  97  97 112  97   1 107   1  97   1  97   1  97   1  97
## [2942]  99  97  46  97  45  97  97  99  97  41  97  99  97  97  97 101  97
## [2959]  97 101  97  97  97  97  97  97  97  97  97  97  97 103  98  97  97
## [2976]  98  98  97  92  97  97  98  97  97  97  97 101  97  92 112 101  41
## [2993]  46 108  46  97  46  97  92  97  98  97  97  97  98  97  97  97  97
## [3010]  97  97  99  98  41  97  67  98  41  97  92  97  46  98   1  98   1
## [3027]  98   1  97   1  97  97  97  97  97  97  92  97  97  98  98  97  97
## [3044]  97  92  97  92  97  97  97  98  97  97  97  97  97  92  97  97  46
## [3061]  97  46  97  41  98  97  97  97  97  97  97  97  97  97  97 102  97
## [3078]  97  97  97  97  97 101  92  97  98  97  97  97  97  97  97 101  97
## [3095]  97  46  97 101  97  97  97  97  99  97  97  97  97 112  98   4  97
## [3112]  98  97  97 112  98  46  98  97  97  97  97  97  97  97  97  41  98
## [3129]   1  97   1  97   1  97   1  97   1  99   9  97   1 101   1  97 101
## [3146]  99  97  97  97  97  97  97  97  98  97  97 101  93  46  97   2 107
## [3163]   1  46   1  97  97  97  46  98   1  97   1  97  98  97  97  97 112
## [3180]  97  92 111  97  97  97  97  97  97  46 100 100  46  97  97  97  97
## [3197]  97 101  97  97  97  97  97  97  97  97  97  97  97  97  97  97  97
## [3214]  97  97  46  97  92  97  97  97 112  97  41  97  46  97  39  97  42
## [3231] 112  98  41  46  92  97  97  46  97  53  46  46  97  44  97  97 101
## [3248]  97  97  97 101  98  97  97 101  97  97  97  98  97  97  97  97  97
## [3265]  97  97 112  97  38  97  36  97  97  97  97  97  97  97  97  97  97
## [3282]  97  97 101  97  97  98  98  92  97  97  98  98  97  97  97  97  97
## [3299]  98  98  97  97 112  98 107 101 112  97 112  97 100  98  97  92  98
## [3316]  97  98  97  97  97  97  97  98  97  97  97 101  97 112  46  41  92
## [3333]  97  92  97  97  97  97  97  97  97  97  97  97  97  97  97  98  97
## [3350] 102 112  97  36  97   1  97  97  97  97  97  98  97  46  98  92  97
## [3367]  97  97  97  92  97  97  97  97  97  97  97  97  46  97  97  97  97
## [3384]  98  97  97  97  97  97  97  97  97  97  97  97  98  46 105  42 101
## [3401]  37  97  98  97  97  97  97  97  98  97  98  97  97  97 111  97  97
## [3418]  97  97  97  97  98  97  97  97  97  97  98  97  97  97  97  97  98
## [3435]  97  98  97  97  97  97  46 101  99  97  97  97  97  97  97  46 112
## [3452]  82  36  46  44   1  46 112  98  42  97  56  98  40  97  46  97  97
## [3469]  97  97  97  98  97  92  98  98  97  47 112 109  34 112 101  40  97
## [3486]  97  97  97  97  97  97  97  97  98  97  97  97  97  97  97  97  97
## [3503]  97  97  97  97  92  97 101  97  97  98  97  97  97  97  99  97  97
## [3520]  97 101  97  97  97  97  97  97  97  97 102  97  97  98  97  97 112
## [3537]  97   4  97   1  97 112  97  40  97   1  97   1  97   1  97  38  97
## [3554]  97  97  92  97  97  97 112  97  39  97  43  97  97  97  97  97  97
## [3571]  97  97  97  97  97  97  97 101 101  92  97  97 101  97  46  97   1
## [3588]  97   1  99   2  97   1  97   1  97   1 112   1  39 112  98  97  92
## [3605]  98  97  97  97 112  99  36  97  40  97  57  97  39  46 104  75  92
## [3622]  46  98  97  46  97 103  97  97  97  97  97  97  97  98 101  97  97
## [3639]  97 110  97  92  97  97  97  97  92  97 112  97  37  97  46  97  44
## [3656]  97  41 101  97  97  97  97  97  97  99  97  97  41 112   7  33   1
## [3673]  97 109  97  41 101 108  97  98  98 112  97  37  97  97  97  97  97
## [3690]  97  98  97  98  97  97  97  97  97  97  97  97  97  97  97  46  97
## [3707]  43  41  44  97  97  97  97  98  98  46  97  98  97  97 102  97  97
## [3724]  99  97  97 112  97  28  97 109  97  97  98  97  97  97  97  97  97
## [3741]  97  97  97  98  97  97  97  97  97  97  97  97  97  97  97  97  97
## [3758]  99  97  97  97  97  92  92  97 101  92 112 101  31  97  46  97  43
## [3775]  97  67  97  41  97  42  97  92  97  97  97  97  97  97 101  97  97
## [3792]  97 112 112  32  50  92 101  46  46  63  97  97  97 112  97  71  97
## [3809]  97  97  97 102  98  97  97  97  98  97  97  97  97  97  97  97  97
## [3826]  97  92  97  97 101  92  97  97  97  97  97  97  97  41  97  57 112
## [3843]  92  37  97  92  98  97  97  97  97  97  97  97  97  97  97  97  97
## [3860]  97  97  97  97  97 101  97  97  97 112  97   2  97   1  97 112  97
## [3877]  46  97  98  97  97 101  97  97  97  97  46  97  41  97  40  97  43
## [3894]  98  97  97  98  97  97 101  98  97  97  97 101  97  97 101 102  97
## [3911]  46  97   2  97   1 112   1  37 105  97  46  97  97  97  97 102  97
## [3928]  97  97 110  97  97  97  98  97  97  97 101  97  97  97 101  97  97
## [3945]  46  97  91 101  41  97  67  97  36  98   1  97 100  97  97  97  97
## [3962]  97  97  97 112  97  46  97  40  46  48  98  46  46  97  99  97  97
## [3979]  97  97  97  97 101  97  97  97  97  46  97 103  97  46  46 101  46
## [3996]  97  41  97  52  97  97  98  97  98  97  92  97  97  97  97  97  97
## [4013] 102  98  97 101 112  97 110 112  97  37  97  92  97  46  97 102  97
## [4030]  97 112  97  97  97  97  98  92  97  98 112  97  36  98  97  97  97
## [4047]  97  97  99  97  97  97  97  97  99  97  97  97  97  97  97  97 101
## [4064]  97  97  97  98  97  98  99  97  97  97  97  98  41  97  93  98  98
## [4081]  97  97  97  97 101  98  46  97  97  97  92  97  97  97  98  97  97
## [4098]  97 110  97  97  97  98  97  97  97  97 101  99  93  97  97  98  97
## [4115]  97  97  97  97  97  97  97  97  97  97  97  97  97  97  41  97   6
## [4132]  97  97  97  97  99  97  97  97  99  97  97  97 100  97  97  97  97
## [4149]  97  97  98  97  99  97  97  97  97  97  97  97 107  97  97  97  97
## [4166]  91  97  97  97  97  97 111  97  97  97 112  98  35  97  98  97  97
## [4183]  97 101  97  46 101   4  97   1  97   1  97 101  46  46  97 112  97
## [4200]   2  97   1 101 112  97  41 102   1  46   1  93  99  97  92  98  97
## [4217]  97  92  97 102  97  97  97  97 112  97  42  98  41  46 111  92  97
## [4234]  97  97  46  97 101  97  98  97  98  98  97  97  97  97  97  97  99
## [4251]  97  97  97  97  98 112  97  36  97  97  97  97  97  97  98  97  97
## [4268]  97  97 102  97  97  97  97  97  97  97  97  97  97  97  97  97  97
## [4285]  98  97  46  98  97  98  97  46  98   1  97   1  41  44  97  97  97
## [4302]  97  97  97  98  97  97  97  97 112  97  19  97 100  97  46  97   1
## [4319]  97   1  97   1  98   2  97  99  97  97  97  97  97  97  97  46  98
## [4336]  97  98  98  46  97 112 103 109  97 112  99  39  97  41  92  97  97
## [4353]  97  46  97 109  97 112 101  39  98 102  97  97  97  92  97  46 110
## [4370]  38  97  41  97 108  97  97  97  97  97  97  97  97  97  98  97  98
## [4387] 112 101 112  97  97  97 101  97  97  97  97  97  97  98  97  98  98
## [4404]  97  98  97  97  98  97  97  97  97  97  97  97  97 103  98  97  97
## [4421]  97  97  97  97 101  92  92  98  97  46  98  97  97  97 112  97  87
## [4438]  97  97  98  98 101  97  97  97 112 101  33  97  99  99  97  98  97
## [4455]  97  97  97  97  98  97  99 112  97   3 112   1   3   1   1   1   1
## [4472] 112   1  36   1   1   1   1   1  97   1  97   1  97   1  97   1  46
## [4489] 112  93  32  97  90  97  97  97  97 107 112  99  41  97  97  97  46
## [4506]  97 101 101 110  97  97  97  97  97  97 112  97  33  98  97  97  97
## [4523]  98  97  97  98  97  97  97  98  97  97  97  97  98  97  98  97  98
## [4540] 112  97  41  98  97  46 103 111  97  97  97 112  92  21  97  63  42
## [4557]  97 101  92  97 101 101  97  97  97  97  97  97  98  97  97  97  97
## [4574]  97  97  97  97  97  97  97  97  97  97  97  46  97  41  97  97  97
## [4591]  97  97 112  97  46  97 101  97  97  97 112  97   2  97   1 101   1
## [4608]  97   1 101   1  97   1  97 112  97  92 100  97  97  97  91  97  98
## [4625]  98  97  46  97  97  97  98  97  97  97  97  97 101  97 101  97  41
## [4642]  97   1  97   1  97   1  97   1  97   2  97   1  97   1  97   1  98
## [4659]   1  93 104  97  97  99  97  97  97  97  97  97  97  97  97  97  97
## [4676]  97  97  98  97  97  97  97  97  97  97  97  97  97  98  97  97  97
## [4693]  92  97  98  42  97  96  97  46  97   1  97   1 110   1  97   1 108
## [4710]   1  97   1  98   3  97   1  97 100  97  46  97   6  97   1  99   1
## [4727]  97  63  97  97  97  97  98  97  98  97  97  97  97  97  97  97  98
## [4744]  98 101  97  97  97  97  97  46  97  42  97  46  97  92  97  98  97
## [4761]  97  97  99  97  97  97  97  97  98 101  97 101  97  97  98  97  97
## [4778]  46  97  61  97  92  97  97  97  97  97 101  46  97  42  97  98  97
## [4795]  97  97  97  97  97 101  97  97  97  98  97  97  99  97  97  97  97
## [4812]  97 105  97  92  97  97 103  97  97  97  97  97  92  97  46  97  92
## [4829]  97  97  97  98  97  97  97  97  97  97  98  97  97 102  97  92  46
## [4846]  97  98  97  97  97  97  97  97  97  97  97  97 102  46  97   5  98
## [4863]  84 101  46  97  47  97  40  98  98  97  97  98  97  46  97  96  98
## [4880]  97  97 101  97 112  98   3  97   1  98 112  97  36  97  97  97  97
## [4897]  97  97 112  97   3  97   1  97   1  97   1  97   1  97   1 101   1
## [4914]  97   1 101   1  97   1  97   1  97  88  97  41  97  97 112  97  37
## [4931]  97  92  97  97  97  97  98  97  97  97  46  97  97  97  98  97 103
## [4948]  97  97  97 103  97  97  98  97  97  98 101  97  97  98  97  97  97
## [4965]  97  97  97  97  97  97  97  97  97  97  97  46  97 102  97  97  97
## [4982]  98  97  97  97  97 112 101   8  97   1  97   1  97   1  98   1  97
## [4999]   1  97
chain <- matrix(c(mu,lambda,phi,beta,delta),nrow=iterations)
colnames(chain) <- c("m","lambda","phi","beta","delta")
plot(mu, type="l", ylab="change point = m")

We see a distinct change between the parameters \(\lambda\) and \(\phi\).

print(mean(mu))
## [1] 87.8142
print(mean(lambda))
## [1] 2.326713
print(mean(phi))
## [1] 0.2211952
plot(phi, type="l", ylab="phi")

plot(lambda, type="l", ylab="lambda")

Now we will examine the distribution of the parameters.

hist(lambda)

hist(phi)

hist(mu)

As we will see from the histograms, \(\lambda\) appears to be a skewed normal distribution, \(\phi\) and \(\mu\) seems to be highly centralized around the mean, with a very high variance.

hist(phi)

hist(lambda)

plot(lambda,phi)

plot(lambda,mu)

plot(beta,delta)

I believe the change point occurs at mean of our \(\mu\) values,87.8142. At this point, see a distinct change in our parameters \(\phi\) and \(\lambda\). The standard deviation of \(\mu\) is: 25.968815 Our 95% confidence interval is: 36, 140. If we plot the number of accidents by year, we can draw a vertical line at the change point and see how the average number of accidents changes. Right of the red line, the number of accidents per year seems to be much lower. From this, we can conclude that our results our consistent with the time-series of the observations.

years <- 1851:1962
change <- 1851+round(mean(mu),0)
plot(years,Y,ylab="accidents")
abline(v=change,col="red")

In Gibbs sampling we are able to sample from the joint distribution (the posterior distribution) by sampling from the full conditional distributions of each parameter. For this, we need to know the posterior distribution. The advantage to Metropolis-Hastings is that we do not need to know the posterior distribution. In this methods, we can choose the proposal distribution, but must be careful to check that the distribution is symmetric and to monitor the acceptance rate in the performance of the algorithm.

mean(mu[b:5000]) + 2*sd(mu[b:5000])
## [1] 139.7545
mean(mu[b:5000]) - 2*sd(mu[b:5000])
## [1] 35.87019

Problem 4

Implement a Metropolis-Hastings simulated annealing optimization to find a least cost solution to the traveling salesman problem (TSP). A sample problem is provided below, which is a 17x17 matrix of cost (C) between every pair of cities in a 17-city TSP problem.

C <- matrix(c(0, 633, 257, 91, 412, 150, 80, 134, 259, 505, 353, 324, 70, 211, 268, 246, 121, 633, 0, 390, 661, 227, 488, 572, 530, 555, 289, 282, 638, 567, 466, 420, 745, 518, 257, 390, 0, 228, 169, 112, 196, 154, 372, 262, 110, 437, 191, 74, 53, 472, 142, 91, 661, 228, 0, 383, 120, 77, 105, 175, 476, 324, 240, 27, 182, 239, 237, 84, 412, 227, 169, 383, 0, 267, 351, 309, 338, 196, 61, 421, 346, 243, 199, 528, 297, 150, 488, 112, 120, 267, 0, 63, 34, 264, 360, 208, 329, 83, 105, 123, 364, 35, 80, 572, 196, 77, 351, 63, 0, 29, 232, 444, 292, 297, 47, 150, 207, 332, 29, 134, 530, 154, 105, 309, 34, 29, 0, 249, 402, 250, 314, 68, 108, 165, 349, 36, 259, 555, 372, 175, 338, 264, 232, 249, 0, 495, 352, 95, 189, 326, 383, 202, 236, 505, 289, 262, 476, 196, 360, 444, 402, 495, 0, 154, 578, 439, 336, 240, 685, 390, 353, 282, 110, 324, 61, 208, 292, 250, 352, 154, 0, 435, 287, 184, 140, 542, 238, 324, 638, 437, 240, 421, 329, 297, 314, 95, 578, 435, 0, 254, 391, 448, 157, 301, 70, 567, 191, 27, 346, 83, 47, 68, 189, 439, 287, 254, 0, 145, 202, 289, 55, 211, 466, 74, 182, 243, 105, 150, 108, 326, 336, 184, 391, 145, 0, 57, 426, 96, 268, 420, 53, 239, 199, 123, 207, 165, 383, 240, 140, 448, 202, 57, 0, 483, 153, 246, 745, 472, 237, 528, 364, 332, 349, 202, 685, 542, 157, 289, 426, 483, 0, 336, 121, 518, 142, 84, 297, 35, 29, 36, 236, 390, 238, 301, 55, 96, 153, 336, 0),ncol=17)

First we will write a function that takes finds the total cost of a solution to the TSP. The function will take in a possible journey - a visit to each of the 17 cities - and the cost matrix. We will test this first by finding the total cost of a journey through all the cities in order, \(1,2,3...17\) and then a randomized journey through all 17.

cities <- 1:17
totalCost <- function(verticies,costMatrix){
     cost <- 0
     for (i in 1:(length(verticies)-1)){
     cost <- cost + costMatrix[verticies[i],verticies[i+1]]
     }
  return(cost)
}
totalCost(cities,C)
## [1] 4601
totalCost(sample(cities),C)
## [1] 3677

Now we will perform simulated annealing that takes three input parameters: the initial temperature \(T0\), \(beta\) (determines annealing schedule), and number of iterations, and return two results: the solution vector \(x\) and the cost of the solution \(s\).

annealing <- function(T0,beta,iterations){
    x <- sample(cities) 
    Sx <- totalCost(x,C)
    xbest <- x 
    sbest <- Sx
    results <- c()

for (i in 1:iterations){
    x <- sample(cities)
    I <- sort(x[1:2])
    Y <- c(x[1:(I[1]-1)],x[(I[2]-1):I[1]],x[I[2]:17])
    Sy = totalCost(Y,C)
    if (Sy < Sx){alpha <- 1} else {alpha <- exp(-(Sy-Sx)/T0)}
    u <- runif(1)
    if (u < alpha){
      x <- Y 
      Sx <- Sy
      }
    T0 <- beta * T0
    xbest <- x
    sbest <- Sx
    results[i] <- sbest
    }
 solutions <- list("vector"=xbest,"cost"=sbest)
 return(solutions)
}
annealing(1,0.9999,10000)
## $vector
##  [1] 15 12 16  3 17  2  7  9  8 14 10 11  1  6  5 13  4
## 
## $cost
## [1] 2845
totalCost(sample(cities),C)
## [1] 4086

As we can see from these two results, the cost of the minimum solution is far less than a random tour through all seven cities. Now we will repeat the cost of a random tour and the simulated annealing, four times each. We see the that the average minimum cost of the annealing process is far less than the average minimum cost of the random tour. The annealing process gives us a different solution vector and minium cost each time. The variance of the minimum cost seems to be quite low, making it a far better option than a random tour.

totalCost(sample(cities),C)
## [1] 4117
totalCost(sample(cities),C)
## [1] 4496
totalCost(sample(cities),C)
## [1] 5003
totalCost(sample(cities),C)
## [1] 4618
annealing(1,0.9999,10000)
## $vector
##  [1] 17  8  5 11 10  4 14 12  3 13  9  7 16  2  6 15  1
## 
## $cost
## [1] 2911
annealing(1,0.9999,10000)
## $vector
##  [1] 17  4 14 15 12 11  9 16  5 13  7  3  8  2  1 10  6
## 
## $cost
## [1] 2655
annealing(1,0.9999,10000)
## $vector
##  [1] 11  2 16  7  3  8 13 10  5  9  1  6 12 15 14 17  4
## 
## $cost
## [1] 2384
annealing(1,0.9999,10000)
## $vector
##  [1] 13  7 14 11 10  3 17  6  4  9  8  5  2 16 12 15  1
## 
## $cost
## [1] 2627

Now we will plot the minimum cost for annealing process against the number of iterations. We will do this on a log scale, for \(10^{i}\) iterations with i=0…5 iterations. We see that the minimum cost decreases greatly each time we raise the power. We have our lowest cost at 100,000 iterations. We can conclude from this example that the performance of the annealing simulations increases greatly with the number of iterations it runs.

power <- 0:5
minCost <- c()
for (i in 0:5){minCost[i+1] <- as.numeric(annealing(1,0.9999,10^i)[[2]])}
plot(power,minCost,xlab="i for 10^i iterations")