The following external R packages are used:
1. Introduction
This is to illustrate changes in marginal probabilities of the “largest value until the current die roll” Markov chain 1 using gganimate package.
2. The largest roll until now
Say you’re rolling a fair six-sided die. Let \(X_0\) be \(0\), and \(X_n\) be the largest value that appears among all of the rolls up to time \(n \geq 1 \text{ }\)1. Clearly, the state space is \(S = \{0, 1, 2, 3, 4, 5, 6 \}\) with initial probabilities \(v_0 = P(X_0 = 0) = 1\), \(v_s = 0\) for all \(s \in S - \{ 0 \}\), and: \[p_{0j} = P(X_{n} = j \text{ | } X_{n - 1} = 0) = \frac{1}{6} \text{ } \forall j \in S - \{ 0\}\]
One fact you can see is that \(X_n\) never decreases by definition, i.e. \(X_{n - 1} \leq X_n\) for \(n \geq 1\). And since we are tossing a fair die, the marginal probability of seeing either side is the same across all sides. So, for example, if \(X_n\) is \(4\), and \(Y_{n + 1}\) is the value that appears at the \((n + 1)\)th roll, then since \(P(Y_{n + 1} = i) = \frac{1}{6}\) for all \(i \in S - \{ 0 \}\), we obtain:
\[\begin{align*} p_{44} &= P(X_{n + 1} = 4 \text{ | } X_n = 4) = P(Y_{n + 1} \leq 4) = \frac{4}{6} \\ p_{45} &= P(X_{n + 1} = 5 \text{ | } X_n = 4) = P(Y_{n + 1} = 5) = \frac{1}{6} \\ &= p_{46} \end{align*}\]
That is, if \(1, 2, 3,\) or \(4\) shows up at the \((n + 1)\)th roll, then the current maximum \(X_{n + 1}\) is still \(4\). Using the same argument for all \(i, j \in S\), we get: \[p_{ij} = \begin{cases} \frac{j}{6} & i = j \\ \frac{1}{6} & i < j \\ 0 & i > j \end{cases}\]
Say \(P_{7 \times 7}\) is the matrix of transition probabilities. Then it should be: \[P_{7 \times 7} = \big[p_{ij} \big]_{i = 0:6, j = 0:6}= \begin{bmatrix} 0 & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} \\ 0 & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} \\ 0 & 0 & \frac{2}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} \\ 0 & 0 & 0 & \frac{3}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} \\ 0 & 0 & 0 & 0 & \frac{4}{6} & \frac{1}{6}& \frac{1}{6} \\ 0 & 0 & 0 & 0 & 0 & \frac{5}{6} & \frac{1}{6} \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}\]
My goal is to visualize the change in marginal probabilities \(P(X_n = i)\), \(i \in S\), as \(n \to \infty\). I will define \(P\)-matrix, and set 50 as the maximum \(n\):
# Settings
max_n <- 50
v <- c(1, rep(0, 6)) # initial probabilities v_0, v_1, ..., v_6
# P matrix
P <- diag(0:6)
for (i in 1:7) {
for (j in 1:7) {
if (i < j) {
P[i, j] <- 1
}
}
}
(P <- P / 6)## [,1] [,2] [,3]
## [1,] 0 0.1666667 0.1666667
## [2,] 0 0.1666667 0.1666667
## [3,] 0 0.0000000 0.3333333
## [4,] 0 0.0000000 0.0000000
## [5,] 0 0.0000000 0.0000000
## [6,] 0 0.0000000 0.0000000
## [7,] 0 0.0000000 0.0000000
## [,4] [,5]
## [1,] 0.1666667 0.1666667
## [2,] 0.1666667 0.1666667
## [3,] 0.1666667 0.1666667
## [4,] 0.5000000 0.1666667
## [5,] 0.0000000 0.6666667
## [6,] 0.0000000 0.0000000
## [7,] 0.0000000 0.0000000
## [,6] [,7]
## [1,] 0.1666667 0.1666667
## [2,] 0.1666667 0.1666667
## [3,] 0.1666667 0.1666667
## [4,] 0.1666667 0.1666667
## [5,] 0.1666667 0.1666667
## [6,] 0.8333333 0.1666667
## [7,] 0.0000000 1.0000000Say \(P^k = P^{(k)}\) stores \(p_{ij}^{(k)} = P(X_{n + k} = j \text{ | } X_n = i)\), i.e. the probability that \(j\) will be the maximum after throwing a die exactly \(k\) times given the current maximum of \(i\). Also, say \(\mu_k\) is the (row) vector of all \(P(X_k = i)\)’s, \(i \in S\). Then given a vector of initial probabilities \(v = \begin{bmatrix} v_0 & v_1 & \dots & v_6 \end{bmatrix}_{1 \times 7}\): \[\mu_k = v P^k\]
Pk <- function(P, k) {
# P: a numeric matrix
# k: a natural number
lst_k <- vector('list', k)
for (i in 1:k) {
lst_k[[i]] <- P
}
Reduce('%*%', lst_k)
}
mu_k <- function(v, P, k) {t(v) %*% Pk(P, k)}Each row of mu_collection consists of \(n\) and \(\mu_n\), in this order, \(n = 0, 1, \dots, 50\):
(mu_collection <-
t(sapply(1:max_n, function(k){mu_k(v, P, k)})) %>%
rbind(v, .) %>%
cbind(0:max_n, .) %>%
'colnames<-'(c('n', 0:6)) %>%
as_tibble())## # A tibble: 51 x 8
## n `0` `1` `2`
## <dbl> <dbl> <dbl> <dbl>
## 1 0 1 0. 0.
## 2 1 0 1.67e-1 1.67e-1
## 3 2 0 2.78e-2 8.33e-2
## 4 3 0 4.63e-3 3.24e-2
## 5 4 0 7.72e-4 1.16e-2
## 6 5 0 1.29e-4 3.99e-3
## 7 6 0 2.14e-5 1.35e-3
## 8 7 0 3.57e-6 4.54e-4
## 9 8 0 5.95e-7 1.52e-4
## 10 9 0 9.92e-8 5.07e-5
## # ... with 41 more rows, and
## # 4 more variables:
## # `3` <dbl>, `4` <dbl>,
## # `5` <dbl>, `6` <dbl>It makes sense to see \(P(X_n = 6) \to 1\) as \(n\) increases. Now we should manipulate mu_collection so that it is ready for plotting a barplot of marginal probabilities at each \(n\):
And here’s how marginal probabilities evolve as \(n \to \infty\):
anim <- ggplot(mu_tidy, aes(x = key, y = value)) +
geom_bar(stat = "identity") +
labs(x = 'Eye', y = 'Probability', title = "n: {closest_state}") +
transition_states(n) +
enter_fade() +
exit_fade()
animate(
anim,
fps = 30, duration = 30,
start_pause = 5,
end_pause = 5
)
Indeed, \(P(X_n = 6) \stackrel{p}{\to} 1\).
Session info
R session info:
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_Canada.1252
## [2] LC_CTYPE=English_Canada.1252
## [3] LC_MONETARY=English_Canada.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_Canada.1252
##
## attached base packages:
## [1] stats graphics
## [3] grDevices utils
## [5] datasets methods
## [7] base
##
## other attached packages:
## [1] tidyr_1.0.0
## [2] gganimate_1.0.4
## [3] ggConvexHull_0.1.0
## [4] dplyr_0.8.3
## [5] reticulate_1.13
## [6] pROC_1.15.3
## [7] ggrepel_0.8.1
## [8] ggplot2_3.2.1
## [9] funpark_0.2.6
## [10] data.table_1.12.6
## [11] boot_1.3-22
## [12] rmarkdown_1.17
## [13] magrittr_1.5
## [14] itertools2_0.1.1
##
## loaded via a namespace (and not attached):
## [1] progress_1.2.2
## [2] prettydoc_0.3.1
## [3] tidyselect_0.2.5
## [4] xfun_0.11
## [5] purrr_0.3.3
## [6] lattice_0.20-38
## [7] colorspace_1.4-1
## [8] vctrs_0.2.0
## [9] htmltools_0.4.0
## [10] yaml_2.2.0
## [11] utf8_1.1.4
## [12] rlang_0.4.2
## [13] pillar_1.4.2
## [14] glue_1.3.1
## [15] withr_2.1.2
## [16] tweenr_1.0.1
## [17] lifecycle_0.1.0
## [18] plyr_1.8.4
## [19] stringr_1.4.0
## [20] munsell_0.5.0
## [21] gtable_0.3.0
## [22] evaluate_0.14
## [23] labeling_0.3
## [24] knitr_1.26
## [25] gifski_0.8.6
## [26] fansi_0.4.0
## [27] Rcpp_1.0.3
## [28] readr_1.3.1
## [29] scales_1.1.0
## [30] backports_1.1.5
## [31] jsonlite_1.6
## [32] farver_2.0.1
## [33] gridExtra_2.3
## [34] png_0.1-7
## [35] hms_0.5.2
## [36] digest_0.6.23
## [37] stringi_1.4.3
## [38] grid_3.6.1
## [39] cli_1.1.0
## [40] tools_3.6.1
## [41] lazyeval_0.2.2
## [42] tibble_2.1.3
## [43] crayon_1.3.4
## [44] pkgconfig_2.0.3
## [45] zeallot_0.1.0
## [46] MASS_7.3-51.4
## [47] ellipsis_0.3.0
## [48] Matrix_1.2-17
## [49] prettyunits_1.0.2
## [50] xml2_1.2.2
## [51] assertthat_0.2.1
## [52] rstudioapi_0.10
## [53] iterators_1.0.12
## [54] R6_2.4.1
## [55] compiler_3.6.1Rosenthal, J. (2019, April 05). STA447/2006 (Stochastic Processes) Lecture Notes, Winter 2019. Retrieved May 23, 2019, from http://probability.ca/jeff/teaching/1819/sta447/notes.pdf↩