将 matlab 循环转换为 R 循环

Question

我有一个转换矩阵如下：

        stage1   stage2 stage3   stage4   stage5
stage1 0.967716 0.017084  0.000 0.000000 0.015200
stage2 0.100000 0.500000  0.200 0.100000 0.100000
stage3 0.200000 0.300000  0.300 0.100000 0.100000
stage4 0.000000 0.000000  0.038 0.917498 0.044502
stage5 0.000000 0.000000  0.000 0.000000 1.000000

下面的matlab代码是这个转移矩阵第一行的条件循环

for i=1:1000
    a=unifrnd(0,1);
    if a<=0.967716 
        stage(i)=1;
    else
         if a<=0.9848
             stage(i)=2;
         else

             stage(i)=5;

         end
    end
end

这意味着在 1000 次迭代中，如果条件为真，均匀分布生成的数字将被分配到五个阶段之一。生成的数字先与停留在当前状态的概率进行比较，如果条件不满足则与上一个概率与转移到下一个状态的概率之和（0.967716 + 0.017084）进行比较，以此类推。

现在我想将这些代码转换为 R 代码。因此，如果我有一个跟踪矩阵“f”，则均匀分布生成的数字将分配给该矩阵的元素之一（列等于上述转换矩阵状态（阶段），行等于迭代）。

f<-matrix(NA, nrow=5, ncol=1000)

...

有人可以帮助我吗？

Answer 1

你的意思是stage(i)=6;而不是stage(6)=6;？

如果是这样，您的 Matlab 代码可以很容易地改写为

stage = discretize(unifrnd(0,1,1000,1),[0,0.92,0.94,0.97,1],[1,2,3,6],'IncludedEdge','right');

R中的代码很相似，你可以在这里找到函数调用 https://www.rdocumentation.org/packages/arules/versions/1.6-4/topics/discretize

Answer 2

如果我很了解您的问题，您希望对马尔可夫过程的多个阶段进行建模。您的过程需要五个阶段。在您的示例中，您只想对概率定律进行建模以从一个阶段传递到另一个阶段。

我接受了您的尝试并对其进行了修改以使其正常工作。

小心：

• 当你想完成矩阵时，你必须同时指定行和列m[i, j]来填充行i和列j的大小写。
• 如果您使用随机性进行模拟，请使用 set.seed 以获得可重复性。
• 在使用像runif 这样的函数之前，请检查它是如何工作的： https://stat.ethz.ch/R-manual/R-devel/library/stats/html/Uniform.html

在这个例子中，马尔可夫过程没有概率停留在同一阶段。

p <- matrix(NA, ncol=5, nrow=50)
set.seed(123) # It fixed the seed for random generation of your computer. Very useful to reproductible code.
for(i in 1:nrow(p)){
  a<-runif(1, 0, 1) # You want only one uniform random number not five
  if(a <= 0.25){
    p[i, 2] <- 2} # p has n rows and you need to allocate one row by one row, so i is needed
  else if(a<= 0.66){
    p[i, 3] <- 3}
  else if(a <= 0.88){
    p[i, 4]<- 4}
  else{
    p[i, 5] <- 5}}

我对您已经尝试的内容进行了一些更改，因为您有五列，您可以在正确的情况下放一列。这个想法是看看你是否找到了概率，这要归功于模拟概率法则来获得每个阶段。就是说：

0.25 to have 2
0.41 to have 3
0.22 to have 4
0.12 to have 5

感谢这段代码：

n <- 10000
p <- matrix(0, ncol=5, nrow=n)
set.seed(123)
for(i in 1:nrow(p)){
  a <- runif(1, 0, 1)
  if(a <= 0.25){
    p[i, 2] <- 1}
  else if(a<= 0.66){
    p[i, 3] <- 1}
  else if(a <= 0.88){
    p[i, 4]<- 1}
  else{
    p[i, 5] <- 1}}

这给出了这些概率：

 colSums(p) / n
 0.0000 0.2475 0.4173 0.2211 0.1141

有效！希望对你有帮助。