Evolution of a binary trait

We simulate the evolution of a binary trait through time

The simulation function

simulate <- function (T, p, lambda) {
  N = 0
  t = 0.0
  state = rbinom(1,1,p)
  states = c(state)
  waitingTimes = c()
  while (t < T) {
    X = rexp(n=1, lambda)
    t = t+X
    if (t < T) {
      N=N+1
      if (state == 0) {
        state = 1
      }
      else {
        state = 0
      }
      states=c(states, state)
      waitingTimes = c(waitingTimes, X)
    }
  }
  return (list(N, states, waitingTimes))
}

We simulate:

Initial values

T = 1.0
p = 0.3
lambda = 10

The actual simulation

elements = simulate(T, p, lambda)
print("Number of substitutions: ")

[1] "Number of substitutions: "

print(elements[[1]])

[1] 12

print("History of the states: ")

[1] "History of the states: "

print(elements[[2]])

 [1] 0 1 0 1 0 1 0 1 0 1 0 1 0

print("Waiting times between substitutions:")

[1] "Waiting times between substitutions:"

print(elements[[3]])

 [1] 0.06138306 0.05248648 0.03947663 0.04999134 0.01849040 0.06142619 0.22715882 0.08331782 0.10947962 0.04423772
[11] 0.03183283 0.03069679

print("Sum of the waiting times:")

[1] "Sum of the waiting times:"

print(sum(elements[[3]]))

[1] 0.8099777