```
= 0.2
lambda
= 0.5
mu
= 10000
n_jumps
= rep(NA, n_jumps + 1)
X_t = rep(NA, n_jumps + 1)
W_n = rep(NA, n_jumps + 1)
T_n
1] = 0
X_t[1] = 0
T_n[
for (n in 2:(n_jumps + 1)){
if (X_t[n - 1] == 0){
- 1] = rexp(1) / lambda
W_n[n = T_n[n - 1] + W_n[n - 1]
T_n[n] = X_t[n - 1] + 1
X_t[n] else {
} - 1] = rexp(1) / (lambda + mu)
W_n[n = T_n[n - 1] + W_n[n - 1]
T_n[n] = X_t[n - 1] + sample(c(-1, 1), 1, prob = c(mu, lambda))
X_t[n]
} }
```

# Birth and Death Chains

## M/M/1 Queue simulation of long run distribution

```
plot(T_n, X_t,
type = "s",
xlab = "t", ylab = "X(t)",
main = "Sample path of M/M/1 queue")
```

```
= data.frame(T_n, X_t, W_n) |>
long_run_distribution slice_head(n = n_jumps) |>
mutate(total_time = sum(W_n)) |>
group_by(X_t) |>
summarize(total_time_in_state = sum(W_n),
total_time = max(total_time),
fraction_time_in_state = total_time_in_state / total_time)
|>
long_run_distribution select(X_t, fraction_time_in_state) |>
kbl(digits = 3) |> kable_styling()
```

X_t | fraction_time_in_state |
---|---|

0 | 0.598 |

1 | 0.240 |

2 | 0.094 |

3 | 0.038 |

4 | 0.017 |

5 | 0.006 |

6 | 0.004 |

7 | 0.001 |

8 | 0.001 |

9 | 0.001 |

10 | 0.000 |

11 | 0.000 |

12 | 0.000 |

13 | 0.000 |

```
|>
long_run_distribution ggplot(aes(x = X_t,
y = fraction_time_in_state)) +
geom_point() +
geom_line() +
labs(x = "State (number in system)",
y = "Long run fraction of time in state")
```

## M/M/1 Queue theoretical stationary distribution

```
= 0.2
lambda
= 0.5
mu
= 0:10
x
= dgeom(x, 1 - lambda / mu)
px
= data.frame(x, px)
stationary_distribution
|>
stationary_distribution kbl(digits = 3) |>
kable_styling()
```

x | px |
---|---|

0 | 0.600 |

1 | 0.240 |

2 | 0.096 |

3 | 0.038 |

4 | 0.015 |

5 | 0.006 |

6 | 0.002 |

7 | 0.001 |

8 | 0.000 |

9 | 0.000 |

10 | 0.000 |

```
|>
stationary_distribution ggplot(aes(x = x,
y = px)) +
geom_point() +
geom_line() +
labs(x = "State (number in system)",
y = "Stationary probability of state")
```

## M/M/1 Queue: Mean time until state 4

```
# Since we want expected time until state 4, treat state 4 as absorbing
# And treat states 0, 1, 2, 3 as transient
# submatrix of transition probability matrix of embedded discrete time chain
# only transition probabilities for "transient" states to "transients" states
= rbind(
QD c(0, 1, 0, 0),
c(5, 0, 2, 0) / 7,
c(0, 5, 0, 2) / 7,
c(0, 0, 5, 0) / 7
)
# Mean amount of time spent in each state before leaving
= c(1 / 0.2, rep(1 / (0.2 + 0.5), 3))
inv_lambda
# Solve system for mean amount of time until state 4
# See end of Handout 23
solve(diag(nrow(QD)) - QD, inv_lambda)
```

`[1] 198.125 193.125 175.625 126.875`

## M/M/1 Queue: Time until state 4

```
= 0.2
lambda
= 0.5
mu
= 4
absorbing_state
= 10000
n_reps
= rep(NA, n_reps)
T_a
for (i in 1:n_reps) {
= 0
x = 0
t_a while(x < absorbing_state) {
= t_a + rexp(1) / (lambda + (x > 0) * mu)
t_a = x + sample(c(-1, 1), 1, prob = c((x > 0) * mu, lambda))
x
}= t_a
T_a[i]
}
hist(T_a)
```

`mean(T_a)`

`[1] 194.5092`