Stepfun function markov

Question

Don't be scared by my long code. What i am wondering is about the last part, the plot(step fun... part. When i enter this into Rstudio i get "stepfun "x" must be ordered increasingly"

Is there any1 here who knows what I have to do to finish this correctly?

bd_process <- function(lambda, mu, initial_state = 0, steps = 100) {
  time_now <- 0
  state_now <- initial_state

  time <- 0
  state <- initial_state
  for (i in 1:steps) { 

if (state_now == 3) {
  lambda_now <- 0
} else {
  lambda_now <- lambda
}

if (state_now == 0) {
  mu_now <- 0
} else {
  mu_now <- mu
}

time_to_transition <- rexp(mu, rate = 1) + rexp(lambda, rate = 1)

X <- rexp(mu, rate = 1)
Y <- rexp(lambda, rate = 1)

if (X < Y) {
  state_now <- state_now - 1
} else {
  state_now <- state_now + 1
}
 time_now <- time_now + time_to_transition 
 time <- c(time, time_now) 
 state <- c(state, state_now) 
}

list(time = time, state = state) }

set.seed(19930628) 

proposal1 <- bd_process(lambda = 2, mu = 10)
proposal2 <- bd_process(lambda = 6, mu = 10)
proposal3 <- bd_process(lambda = 10, mu = 10)


time1 <- proposal1$time
state1 <- proposal1$state

plot(stepfun(time1[-1], state1), 
 do.points = FALSE,
 xlab = "Tid",
 ylab = "Tillstånd",
 main = "",
 yaxt = "n")
axis(2, at = c(0, 1, 2, 3), las = 2)

Well, the error message sounds like the `stepfun` function expects your data to be sorted. Have you tried that? Maybe read the help page for `?stpfun`. It's a bit difficult to help when you don't describe what the desired outcome is at all. No idea what "correct" means to you. — MrFlick, Jul 17 '16 at 18:07
The question is really long & it got voted down the last time i tried to explain it. — PeterNiklas, Jul 17 '16 at 18:25
I would just want to know what i need to change, so that i can KNIT. For instance I´ve tried time1[-1], state1 and time1[100], state1[100] and so on. — PeterNiklas, Jul 17 '16 at 18:26
If you have transitions at random intervals, then you need to cumsum the interval length to create a time variable for plotting. — IRTFM, Jul 17 '16 at 21:46
You also need to figure out why: `length(time1) == 1001` and `length(state1) == 101` — IRTFM, Jul 17 '16 at 22:24

IRTFM · Accepted Answer · 2016-07-18T15:40:46.257

I don't know what your code is doing but you've asked us not to worry about that. At the moment it appears that you have only constructed "time intervals" but now need to "stack them together" or "integrate" them along a proper time axis. In order to plot a simulation of a stepfunction, you should be using cumsum to construct an increasing time1 vector. Because the "time" and "state" variables are of such different lengths a quick fix to the function arguments is trimming the time1 vector so it is the correct length for the state1 variable, and you get no error with:

plot(stepfun(cumsum(time1[2:101]), state1), 
 do.points = FALSE,
 xlab = "Tid",
 ylab = "Tillstånd",
 main = "",
 yaxt = "n")
axis(2, at = c(0, 1, 2, 3), las = 2)

Maybe if you "march step-by-step" through the code and explain the code (to yourself and the rest of us) using comments you will figure out why you have 10 times as many time1's as you have state1's. I suspect it may have something to do with using "mu" as the first argument to rexp(mu, rate = 1). The first argument to random number generators in R is usually a positive integer that determines length (the number of random numbers) from the distribution.

i know this might be hard to answer with limited information, but is it possible to define the time1 vektorn with ascending order? — PeterNiklas, Jul 20 '16 at 17:11
If the "times" are not successive and there is no censoring or domination of intermediate times by any later "times" then there is both a `sort` and an `order` function. Generally the `order` function would be preferred since it is more likely in simulations that multiple vectors would need to have the the same re-ordering operation. If you just `sort` one vector and not the others, then any temporal association will be lost. — IRTFM, Jul 20 '16 at 17:48

Stepfun function markov

1 Answers1