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

Question

У меня есть следующий код MATLAB, и я работаю над переводом его на 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);

Это Пуассоновский процесс с точки зрения процессов обновления. Я сделал все возможное в R, но что-то не так в моем коде:

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

Цикл печатает произвольное количество массивов, и только последний сохраняется на tarr после него.

Результат должен выглядеть так ...

Заранее спасибо. Все интересные и поддерживающие комментарии будут вознаграждены.

Answer 1

В дополнение к предыдущему комментарию в сценарии 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 и не подходит для перезаписи!

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

Answer 2

Сегодня я скачал 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 = "")

Answer 3

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

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

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

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

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

FALSE == 1 # возвращает FALSE

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

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

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

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

vector <- c(vector,newvector)

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

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

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

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

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

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

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