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 代码，看看发生了什么。

前 4 行为变量赋值。
第五行创建一个包含 40 (nproc) 个零的数组。
第六行（后面好像没用到）创建一个模式为'list'的空向量。
第七行开始一个while循环。我怀疑这条线应该说 而 tarr 的最小值小于或等于 T ... 或者它应该说 while i 小于或等于 T ... 它实际上取单个布尔值的最小值 (tarr[i] <= T)。现在这可以工作了，因为 TRUE 和 FALSE 被视为数字。即：

TRUE == 1 # returns TRUE

FALSE == 0 # returns TRUE

TRUE == 0 # returns FALSE

FALSE == 1 # returns FALSE

然而，由于 tarr[i] 的值取决于一个随机数（见第 8 行），这可能导致相同的代码运行ning 每次执行时都不同。这或许可以解释为什么代码“打印出任意数量的数组”。

第八行似乎用右边的计算覆盖了tarr的赋值。因此它采用 tarr[i] 的单个值并从中减去运行if(proc) 除以 4 (lambda) 的自然对数——这给出了 40 个不同的值。这 40 个来自上次循环的不同值存储在 tarr 中。如果你想在循环中每次都存储所有四十个值，我建议将它存储在矩阵或数据框中。如果那是你想要做的，这里有一个将它存储在矩阵中的例子：

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

有关详细信息，请参阅答案。此外，由于它是每个运行的一组值，您可以将它们保存在一个向量中并简单地附加到每个循环的向量中，如下所示：

vector <- c(vector,newvector)

第九行将i加1
第十行打印tarr.
第十一行结束循环语句。
然后循环后tarr2赋值1/tarr。同样，这将是上次循环中的 40 个值（第 8 行）
然后X赋给tarr2的最小值。
最后一行绘制了这个单一值。

另请注意 runif 来自均匀分布的样本——如果您正在寻找泊松分布，请参阅：Poisson

希望对您有所帮助！让我知道是否还有更多我可以帮助的事情。

Answer 2

添加到之前的评论中，在 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() 将 return 矩阵的最小值（只有一个数字）和 M: min() return s 每列的最小值。因此，对于 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 3

我今天下载了 GNU Octave 来实际运行 MatLab 代码。在查看代码运行ning 之后，我花了几个星期的时间来回答 @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 = "")

R中的泊松过程算法（更新过程的角度）

Poisson Process algorithm in R (renewal processes perspective)

r

loops

simulation

stochastic