Loop in a dataset simulation

Question

I hope to get help on the following problem in R.

I have the folowing code to generate 30 column dataset based on an exponential distribuition:

x0=0
xmax=8000
xout=3000
lambda=0.0002
n=1

  x1=x0+rexp(n,lambda)-xout
  x2=x1+rexp(n,lambda)-xout
  x3=x2+rexp(n,lambda)-xout
  x4=x3+rexp(n,lambda)-xout
  x5=x4+rexp(n,lambda)-xout
  x6=x5+rexp(n,lambda)-xout
  x7=x6+rexp(n,lambda)-xout
  x8=x7+rexp(n,lambda)-xout
  x9=x8+rexp(n,lambda)-xout
  x10=x9+rexp(n,lambda)-xout
  x11=x10+rexp(n,lambda)-xout
  x12=x11+rexp(n,lambda)-xout
  x13=x12+rexp(n,lambda)-xout
  x14=x13+rexp(n,lambda)-xout
  x15=x14+rexp(n,lambda)-xout
  x16=x15+rexp(n,lambda)-xout
  x17=x16+rexp(n,lambda)-xout
  x18=x17+rexp(n,lambda)-xout
  x19=x18+rexp(n,lambda)-xout
  x20=x19+rexp(n,lambda)-xout
  x21=x20+rexp(n,lambda)-xout
  x22=x21+rexp(n,lambda)-xout
  x23=x22+rexp(n,lambda)-xout
  x24=x23+rexp(n,lambda)-xout
  x25=x24+rexp(n,lambda)-xout
  x26=x25+rexp(n,lambda)-xout
  x27=x26+rexp(n,lambda)-xout
  x28=x27+rexp(n,lambda)-xout
  x29=x28+rexp(n,lambda)-xout
  x30=x29+rexp(n,lambda)-xout

I have three doubts:

1 - Is there any way to write this function in a reduced form?

2 - This row (30 columns) needs to be simulated 10,000 times. How to do this in a loop?

3 - The values of each cell (x1, x2, x3 ...) must be limited to the interval x0 and xmax (0-8000). How to do this?

score 2 · Answer 1 · answered May 13 '19 at 17:48

That depends on what you want to do with values over 8000. Here's a solution that just takes those values and wraps them around with a modulo operator.

library(tidyverse)

test <- data.frame(x0 = rep(0, n))
for (i in 1:30) {
  new_col <- sym(paste0("x", i))
  old_col <- sym(paste0("x", i - 1))
  test <- test %>% 
    mutate(!!new_col := (!!old_col + rexp(n, lambda) - xout) %% xmax)
}

I don't know how familiar you may or may not be with the tidyverse and tidy evaluation, which I've used here liberally. The !! operator, combined with sym(), turns the variable names into actual variables. The %>% operator "pipes" data from one function to the next. The := operator is needed only if you want to make assignments with a !! on the lefthand side.

I think this is my first time actually trying to post an answer on StackOverflow, so be easy on me! :)

Thank you @Gregor! I see you're well-past the first-answer mark. ;) — Benjamin, May 13 '19 at 18:00

Gregor Thomas · Answer 2 · 2019-05-13T17:51:57.003

Is there any way to write this function in a reduced form?

I would do it like this. Pretty sure this is equivalent.

ncol = 30

row = rexp(ncol, lambda)
row = cumsum(row) - xout * (1:ncol)

This row (30 columns) needs to be simulated 10,000 times. How to do this in a loop?

Use replicate with the code above:

sim_data = t(replicate(10000, {
  row = rexp(ncol, lambda)
  row = cumsum(row) - xout * (1:ncol)
}))

replicate gives 10000 columns and 30 rows. We use t() to transpose it to 10000 rows with 30 columns.

The values of each cell (x1, x2, x3 ...) must be limited to the interval x0 and xmax (0-8000). How to do this?

Use pmin() and pmax(). Not sure if you want this done before or after the cumulative summing...

sim_data = t(replicate(10000, {
  row = rexp(ncol, lambda)
  row = cumsum(row) - xout * (1:ncol)
  row = pmax(0, row)
  row = pmin(xmax, row)
  row
}))

score 1 · Accepted Answer · answered May 13 '19 at 18:20

As I'm fairly new to R myself, I thought it would be good practice to try to write this out. Perhaps not the most efficient code, but it works:

xmax <- 8000
xout <- 3000
lambda <- 0.0002
n <- 1
iterations <- 30

df <- data.frame(matrix(ncol = 31, nrow = iterations))
names(df) <- c(paste("x", 0:30, sep=""))

for (j in 1:iterations) {
  df$x0[j] <- 0
  df$x1[j] <- df$x0[j] + rexp(n,lambda)-xout
  if (df$x1[j] < 0) {
    df$x1[j] <- 0
  }
  if (df$x1[j] > 8000) {
    df$x1[j] <- 8000
  }
  for (i in 3:31) {
    df[j,i] <- df[j, i-1] + rexp(n,lambda)-xout
    if (df[j,i] < 0) {
      df[j,i] <- 0
    }
    if (df[j,i] > 8000) {
      df[j,i] <- 8000
    }
  }
}

You can change iterations to 30000, for testing purposes I've used 30. Also I didn't know if you wanted to limit to 0 and 8000 before or after the next iterations, I've done it before.

Loop in a dataset simulation

3 Answers3