inst/doc/stat2003-vignette.R

In paulnorthrop/stat2003: stat2003: Introduction to Applied Probability at UCL

## ----set-options, echo=FALSE, cache=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
options(width = 999)

## ---- include = FALSE----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
knitr::opts_chunk$set(comment = "#>", collapse = TRUE)

## ----setup, include=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
devtools::load_all() # reload all code (after saving them)   or Ctrl-shift-L

## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
library(stat2003)
P <- matrix(c(0, 1, 0, 0, 0, 0, 0,
              0, 0, 1, 0, 0, 0, 0,
              0, 1/3, 0, 2/3, 0, 0, 0,
              0, 0, 0, 0, 2/3, 1/3, 0,
              0, 2/3, 0, 1/3, 0, 0, 0,
              0, 0, 0, 0, 0, 1/3, 2/3,
              0, 0, 0, 0, 0, 1/4, 3/4), nr = 7, nc=7, byrow = T)

# input the transition matrix, initial distribution and statespace to
# create a new discrete-time finite statespace Markov chain in the class stat2003.d
MC<- new("stat2003.d", p_start = c(1/7, 1/7, 1/7, 1/7, 1/7, 1/7, 1/7), p = P,
         statespace = LETTERS[1 : 7] )


## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
dmc_irreclass(MC)

## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
dmc_inv(MC, show_mat = FALSE)
dmc_equi(MC)

## ---- fig.align='center'-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
dmc_tp(MC, 10, graph = TRUE)

## ---- fig.align='center'-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
dmc_simu(MC, 30)

## ---- fig.align='center'-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
simu <- dmc_simu(MC, 10000)
table(simu[5000:10000]) / sum(table(simu[5000:10000]))

## ---- fig.align='center'-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# we have to input this Markov chain into the class stat2003.d first.
p <- matrix(c(1/3, 0, 2/3, 0,
             0, 0, 0.5, 0.5,
             1, 0, 0, 0,
             0, 1/4, 0, 3/4), nr = 4, byrow = TRUE)

# Note that, p_start is set to (0, 1, 0, 0), which is corresponding to P(X_0 = 2) = 1 in part 4 of the question.
A <- new("stat2003.d", p_start = c(0, 1, 0, 0), p = p,
          statespace = as.character(1:4))

#(a)
dmc_irreclass(A)

#(b)
dmc_inv(A, show_mat = FALSE)

#(c)
# From the lecture slides
# we know that the mean recurrence time of state 1 is 
# the reciprocal of the first element of invariant distribution
# if the Markov chain is an irreducible ergodic Markov chain.

# from the answer of (a), we know that 
# this Markov chain is an irreducible ergodic Markov chain,
# also based on the answer of (b),
# mean recurrence time of state 1 is
1/0.6

#(d)
# we have already set the initial distribution to (0, 1, 0, 0),
# therefore, we can calculate it directly
dmc_tp(A, 4)
# answer is 0.0859375 

## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
BD1 <- new("stat2003.bd", lb = 2.5, ld = 3, im = 5, em = 1,
            nstate = 6)

## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
bd_rate(BD1)
bd_trans(BD1)

## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
bd_equi(BD1)

## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
BD2 <- new("stat2003.bd", lb = 2.5, ld = 3, im = 0, em = 1,
            nstate = 6)
bd_rate(BD2)
bd_trans(BD2)
bd_equi(BD2)

## ---- fig.align='center'-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
bd_simu(BD1, 15, ini_state = 4, label = 5)
bd_simu(BD2, 15, ini_state = 4, label = 5)

## ---- fig.align='center'-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
BD3 <- new("stat2003.bd", lb = 2.5, ld = 3, im = 5, em = 1,
            nstate = 0)
simu <- bd_simu(BD3, 15, ini_state = 4, label = 5)