Generate Markov Trajectories

Simulate individuals from a Markov process defined by a transition matrix, time spent in each time and initial probabilities.

generate_Markov(
  n = 5,
  K = 2,
  P = (1 - diag(K))/(K - 1),
  lambda = rep(1, K),
  pi0 = c(1, rep(0, K - 1)),
  Tmax = 1,
  labels = NULL
)

Arguments

n

number of trajectories to generate

K

number of states

P

matrix containing the transition probabilities from one state to another. Each row contains positive reals summing to 1.

lambda

time spent in each state

pi0

initial distribution of states

Tmax

maximal duration of trajectories

labels

state names. If NULL, integers are used

Value

a data.frame with 3 columns: id, id of the trajectory, time, time at which a change occurs and state, new state.

Details

For one individual, assuming the current state is \(s_j\) at time \(t_j\), the next state and time is simulated as follows:

1. generate one sample, \(d\), of an exponential law of parameter lambda[s_j]

2. define the next time values as: \(t_{j+1} = t_j + d\)

3. generate the new state \(s_{j+1}\) using a multinomial law with probabilities Q[s_j,]

Author

Cristian Preda

Examples

# Simulate the Jukes-Cantor model of nucleotide replacement
K <- 4
PJK <- matrix(1 / 3, nrow = K, ncol = K) - diag(rep(1 / 3, K))
lambda_PJK <- c(1, 1, 1, 1)
d_JK <- generate_Markov(
  n = 100, K = K, P = PJK, lambda = lambda_PJK, Tmax = 10,
  labels = c("A", "C", "G", "T")
)

head(d_JK)
#>   id      time state
#> 1  1 0.0000000     A
#> 2  1 0.2092139     T
#> 3  1 1.7666995     A
#> 4  1 2.1512726     C
#> 5  1 2.6233965     A
#> 6  1 3.2061763     G