Poisson Process algorithm in R (renewal processes perspective)

Question

I have the following MATLAB code and I'm working to translating it to R:

nproc=40
T=3
lambda=4
tarr = zeros(1, nproc);
i = 1;
while (min(tarr(i,:))<= T)
tarr = [tarr; tarr(i, :)-log(rand(1, nproc))/lambda];
i = i+1;
end
tarr2=tarr';
X=min(tarr2);
stairs(X, 0:size(tarr, 1)-1);

It is the Poisson Process from the renewal processes perspective. I've done my best in R but something is wrong in my code:

nproc<-40
T<-3
lambda<-4
i<-1
tarr=array(0,nproc)
lst<-vector('list', 1)
while(min(tarr[i]<=T)){
tarr<-tarr[i]-log((runif(nproc))/lambda)
i=i+1
print(tarr)
}
tarr2=tarr^-1
X=min(tarr2)
plot(X, type="s")

The loop prints an aleatory number of arrays and only the last is saved by tarr after it.

The result has to look like...

Thank you in advance. All interesting and supportive comments will be rewarded.

Answer 1

I downloaded GNU Octave today to actually run the MatLab code. After looking at the code running, I made a few tweeks to the great answer by @Croote

nproc <- 40
T0 <- 3
lambda <- 4
i <- 1

tarr <- matrix(rep(0, nproc), nrow = 1, ncol = nproc)

while(min(tarr[i, ]) <= T0){

  temp <- matrix(tarr[i, ] - log(runif(nproc))/lambda, nrow = 1) #fixed paren
  tarr <- rbind(tarr, temp)
  i = i + 1
}

tarr2 = t(tarr) #takes transpose

library(ggplot2)
library(Rfast)
min_for_each_col <- colMins(tarr2, value = TRUE)
qplot(seq_along(min_for_each_col), sort(min_for_each_col), geom="step")

Edit: Some extra plotting tweeks -- seems to be closer to the original

qplot(seq_along(min_for_each_col), c(1:length(min_for_each_col)), geom="step", ylab="", xlab="")
#or with ggplot2
df1 <- cbind(min_for_each_col, 1:length(min_for_each_col)) %>% as.data.frame
colnames(df1)[2] <- "index"
ggplot() +
  geom_step(data = df1, mapping = aes(x = min_for_each_col, y = index), color = "blue") +
    labs(x = "", y = "")

Answer 2

Adding on to the previous comment, there are a few things which are happening in the matlab script that are not in the R:

• [tarr; tarr(i, :)-log(rand(1, nproc))/lambda]; from my understanding, you are adding another row to your matrix and populating it with tarr(i, :)-log(rand(1, nproc))/lambda]. You will need to use a different method as Matlab and R handle this type of thing differently.
• One glaring thing that stands out to me, is that you seem to be using R: tarr[i] and M: tarr(i, :) as equals where these are very different, as what I think you are trying to achieve is all the columns in a given row i so in R that would look like tarr[i, ]
• Now the use of min is also different as R: min() will return the minimum of the matrix (just one number) and M: min() returns the minimum value of each column. So for this in R you can use the Rfast package Rfast::colMins.
• The stairs part is something I am not familiar with much but something like ggplot2::qplot(..., geom = "step") may work.

Now I have tried to create something that works in R but am not sure really what the required output is. But nevertheless, hopefully some of the basics can help you get it done on your side. Below is a quick try to achieve something!

nproc <- 40
T0 <- 3
lambda <- 4
i <- 1

tarr <- matrix(rep(0, nproc), nrow = 1, ncol = nproc)

while(min(tarr[i, ]) <= T0){
# Major alteration, create a temporary row from previous row in tarr
  temp <- matrix(tarr[i, ] - log((runif(nproc))/lambda), nrow = 1)
# Join temp row to tarr matrix
  tarr <- rbind(tarr, temp)
  i = i + 1
}
# I am not sure what was meant by tarr' in the matlab script I took it as inverse of tarr
# which in matlab is tarr.^(-1)??
tarr2 = tarr^(-1)

library(ggplot2)
library(Rfast)
min_for_each_col <- colMins(tarr2, value = TRUE)
qplot(seq_along(min_for_each_col), sort(min_for_each_col), geom="step")

As you can see I have sorted the min_for_each_col so that the plot is actually a stair plot and not some random stepwise plot. I think there is a problem since from the Matlab code 0:size(tarr2, 1)-1 gives the number of rows less 1 but I cant figure out why if grabbing colMins (and there are 40 columns) we would create around 20 steps. But I might be completely misunderstanding! Also I have change T to T0 since in R T exists as TRUE and is not good to overwrite!

Hope this helps!

Answer 3

I'm not too familiar with renewal processes or matlab so bear with me if I misunderstood the intention of your code. That said, let's break down your R code step by step and see what is happening.

• The first 4 lines assign numbers to variables.
• The fifth line creates an array with 40 (nproc) zeros.
• The sixth line (which doesnt seem to be used later) creates an empty vector with mode 'list'.
• The seventh line starts a while loop. I suspect this line is supposed to say while the min value of tarr is less than or equal to T ... or it's supposed to say while i is less than or equal to T ... It actually takes the minimum of a single boolean value (tarr[i] <= T). Now this can work because TRUE and FALSE are treated like numbers. Namely:

TRUE == 1 # returns TRUE

FALSE == 0 # returns TRUE

TRUE == 0 # returns FALSE

FALSE == 1 # returns FALSE

However, since the value of tarr[i] depends on a random number (see line 8), this could lead to the same code running differently each time it is executed. This might explain why the code "prints an aleatory number of arrays ".

• The eight line seems to overwrite the assignment of tarr with the computation on the right. Thus it takes the single value of tarr[i] and subtracts from it the natural log of runif(proc) divided by 4 (lambda) -- which gives 40 different values. These fourty different values from the last time through the loop are stored in tarr. If you want to store all fourty values from each time through the loop, I'd suggest storing it in say a matrix or dataframe instead. If that's what you want to do, here's an example of storing it in a matrix:

for(i in 1:nrow(yourMatrix)){
//computations
yourMatrix[i,] <- rowCreatedByComputations
}

See this answer for more info about that. Also, since it's a set number of values per run, you could keep them in a vector and simply append to the vector each loop like this:

vector <- c(vector,newvector)

• The ninth line increases i by one.

• The tenth line prints tarr.

• the eleveth line closes the loop statement.

• Then after the loop tarr2 is assigned 1/tarr. Again this will be 40 values from the last time through the loop (line 8)

• Then X is assigned the min value of tarr2.

• This single value is plotted in the last line.

Also note that runif samples from the uniform distribution -- if you're looking for a Poisson distribution see: Poisson

Hope this helped! Let me know if there's more I can do to help.