Алгоритм процесса Пуассона в R (перспектива процессов обновления)

Сегодня я скачал GNU Octave, чтобы запустить код MatLab. Посмотрев на работающий код, я внес несколько изменений в отличный ответ @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")

Редактировать: некоторые дополнительные настройки сюжета - кажется, ближе к оригиналу

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 = "")

В дополнение к предыдущему комментарию в сценарии Matlab происходит несколько вещей, которых нет в R:

• [tarr; tarr(i, :)-log(rand(1, nproc))/lambda]; Насколько я понимаю, вы добавляете еще одну строку в свою матрицу и заполняете ее tarr(i, :)-log(rand(1, nproc))/lambda]. Вам нужно будет использовать другой метод, поскольку Matlab и R обрабатывают этот тип вещей по-разному.
• Одна бросающаяся в глаза вещь, которая выделяется для меня, заключается в том, что вы, кажется, используете R: tarr[i] и M: tarr(i, :) как равные, где они очень разные, поскольку я думаю, что вы пытаетесь достичь, это все столбцы в данной строке i, поэтому в R это было бы похож на tarr[i, ]
• Теперь использование min также отличается, поскольку R: min() возвращает минимум матрицы (всего одно число), а M: min() возвращает минимальное значение каждого столбца. Так что для этого в R можно использовать пакет Rfast Rfast::colMins.
• Часть stairs - это то, с чем я мало знаком, но что-то вроде ggplot2::qplot(..., geom = "step") может сработать.

Теперь я попытался создать что-то, что работает в R, но я не уверен, что на самом деле требуется. Но, тем не менее, мы надеемся, что некоторые из основ помогут вам сделать это на вашей стороне. Ниже приведена быстрая попытка чего-то добиться!

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")

Как видите, я отсортировал min_for_each_col так, чтобы график был на самом деле ступенчатым, а не каким-то случайным ступенчатым графиком. Я думаю, что есть проблема, так как из кода Matlab 0:size(tarr2, 1)-1 дает количество строк меньше 1, но я не могу понять, почему, если захватить colMins (а есть 40 столбцов), мы создадим около 20 шагов. Но я могу быть совершенно не понимаю! Также я изменил T на T0, так как в R T существует как TRUE и его нельзя перезаписывать!

Надеюсь это поможет!

Я не слишком хорошо знаком с процессами обновления или Matlab, так что потерпите меня, если я неправильно понял намерение вашего кода. Тем не менее, давайте шаг за шагом разберем ваш код R и посмотрим, что происходит.

• Первые 4 строки присваивают номера переменным.
• Пятая строка создает массив из 40 (nproc) нулей.
• Шестая строка (которая, похоже, не будет использоваться позже) создает пустой вектор с режимом «список».
• Седьмая строка запускает цикл while. Я подозреваю, что в этой строке должно быть указано в то время как минимальное значение tarr меньше или равно T ... или должно быть указано в то время как i меньше или равно T .. . На самом деле он принимает минимум одного логического значения (tarr[i] <= T). Теперь это может работать, потому что TRUE и FALSE обрабатываются как числа. А именно:

ИСТИНА == 1 # возвращает ИСТИНА

ЛОЖЬ == 0 # возвращает ИСТИНА

ИСТИНА == 0 # возвращает ЛОЖЬ

ЛОЖЬ == 1 # возвращает ЛОЖЬ

Однако, поскольку значение tarr[i] зависит от случайного числа (см. строку 8), это может привести к тому, что один и тот же код будет работать по-разному при каждом его выполнении. Это может объяснить, почему код выводит произвольное количество массивов.

• Восьмая строка, кажется, перезаписывает назначение tarr вычислением справа. Таким образом, он берет единственное значение tarr[i] и вычитает из него натуральный логарифм runif(proc), деленный на 4 (лямбда), что дает 40 различных значений. Эти сорок различных значений, полученных в последний раз в цикле, хранятся в tarr. Если вы хотите сохранить все сорок значений каждый раз в цикле, я бы предложил вместо этого сохранить их, скажем, в матрице или фрейме данных. Если это то, что вы хотите сделать, вот пример хранения в матрице:

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

См. этот ответ для получения дополнительной информации об этом. Кроме того, поскольку это заданное количество значений за прогон, вы можете хранить их в векторе и просто добавлять к вектору каждый цикл следующим образом:

vector <- c(vector,newvector)

• Девятая строка увеличивает i на единицу.

• Десятая строка выводит tarr.

• одиннадцатая строка закрывает оператор цикла.

• Затем после цикла tarr2 присваивается 1/tarr. Опять же, это будет 40 значений, полученных в последний раз в цикле (строка 8).

• Затем X присваивается минимальное значение tarr2.

• Это единственное значение отображается в последней строке.

Также обратите внимание, что runif образцы из единого дистрибутива — если вы ищете дистрибутив Пуассона, см.: Пуассон

Надеюсь, это помогло! Дайте мне знать, если я могу чем-то помочь.