从Matlab代码转换为R代码的矩阵

Question

这段 Matlab 代码是为了根据生态系统中物种损失来计算生态系统功能概率而创建的。现在，这段代码需要被翻译成 R 语言。但我在将 Matlab 中的矩阵操作转换为 R 时遇到了问题。

在 Matlab 中，以下是我尝试翻译成 R 代码的部分：

for j=1:N+1
multi_matrix4(:,j)=matrix(:,1);
end

在R中，我将这段代码放在了for循环里：

+ multi.matrix4 <- matrix[,1,drop=FALSE]
+ multi.matrix4 <- multi.matrix4[,j,drop=FALSE]
+ class(multi.matrix4)

这是来自R的消息，位于for循环下面：

Error: subscript out of bounds

我的问题是： 如何使用R进行这种矩阵操作？

没有最后一个图形的Matlab代码如下：

clear all

% No of permutations
sim=1000;

% Total No of ecosystem functions    
N=3;

%Total dimensions
J=3;

% Total No of species in pool
total_species=4;

% No of species drawn from pool
species=4;
multi_matrix=zeros(total_species,N);

% "Threshold"
t=.5;

result=zeros(sim,J);

for i=1:sim

% %Uniformly increasing trait values
for j=1:N
matrix=rand(total_species,2);
matrix(:,1)=linspace(0,1,total_species);
matrix=sortrows(matrix,2);
multi_matrix4(:,j)=matrix(:,1);
end

%Complete covariance
matrix=rand(total_species,2);
matrix(:,1)=linspace(0,1,total_species);
matrix=sortrows(matrix,2);
for j=1:N+1
multi_matrix4(:,j)=matrix(:,1);
end

% Excess of high trait values
for j=1:N
matrix=rand(total_species,2);
X=1:total_species;X=X';
matrix(:,1)=1-exp(-0.02*X.^2);
matrix=sortrows(matrix,2);
multi_matrix4(:,j)=matrix(:,1);
end


% Deficiency of high trait values
for j=1:N
matrix=rand(total_species,2);
X=1:total_species;X=X';
% matrix(:,1)=exp((X./22.6).^3)-1;
matrix(:,1)=exp((X./13.55).^3)-1;
matrix=sortrows(matrix,2);
multi_matrix4(:,j)=matrix(:,1);
end


% Reading empirical data
warning off
% [NUMERIC,txt]=xlsread('Plant_6.xls','Sheet1');
Exp07_2 = [ 0 0.72 0.70 ; 1 1 0 ; 0.62 0 1 ; 0.36 0.69 0.61]
multi_matrix(1:total_species,1:N)=Exp07_2;
random=rand(1,N);
multi_matrix(total_species+1,1:N)=random;
multi_matrix2=sortrows(multi_matrix',total_species+1);
multi_matrix3=multi_matrix2';
multi_matrix4=multi_matrix3(1:total_species,:);
warning on


    % adding a sorting column
    random2=rand(total_species,1);
    multi_matrix4(:,N+1)=random2;
    sort_multi_matrix=sortrows(multi_matrix4,N+1);

    % loop adding one function at a time
    for j=1:J

        loss_matrix=sort_multi_matrix(1:species,1:j);
        max_value=loss_matrix>=t;
        B=any(max_value',2);
        C=all(B);
        result(i,j)=sum(C);

    end

end

% reporting
res=mean(result);
res'

这段 R 代码如下：

rm()

#No of permutation
sims <- 1000;

#Total number of ecosystem functions
N <- 3

#Total dimensions
J <- 3

#Total number of species in pool
total.species <- 4

#No of species drawn from pool
species <- 4

multi.matrix <- matrix(0, nrow=total.species, ncol=N)
class(multi.matrix)

# $Threshold$
t <- .5;

# The results are to be put in a matrix
result <- matrix(0, nrow=sims, ncol=J)

for (i in 1 : sims)
{

#Uniformly increasing trait values
for (j in 1 : N)
{
matrix <- matrix(runif(total.species*2),total.species)
class(matrix)
matrix[,1] <- seq(0,1, len=total.species) # test 2
class(matrix)
matrix <- matrix[order(matrix( ,2)),]
class(matrix)
# multi.matrix4[,j,drop=FALSE] = matrix[,1,drop=FALSE]
multi.matrix4 <- matrix[,1,drop=FALSE]
multi.matrix4 <- multi.matrix4[,j,drop=FALSE]
class(multi.matrix4)
}

# Complete covariance
matrix <- matrix(runif(total.species*2),total.species)
class(matrix)
matrix[,1] <- seq(0, 1, len=total.species)
class(matrix)
matrix <- matrix[order(matrix( ,2)),]
class(matrix)
for (j in 1 : N + 1)
{multi.matrix4 <- matrix[,1,drop=FALSE]
multi.matrix4 <- multi.matrix4[,j,drop=FALSE]
class(multi.matrix4)
}

# Excess of high trait values
for (j in 1 : N)
{matrix <- matrix(runif(total.species*2),total.species)
class(matrix)
X <- 1 : total.species
X <- t(X)
matrix[,1] <- c(1 - exp(-0.02 %*% X^2)) # Hie... p. 8
matrix <- matrix[order(matrix( ,2)),]
# multi.matrix4[,j,drop=FALSE] <- matrix[,1,drop=FALSE]
# multi.matrix4[,j,drop=FALSE] <- matrix[,1]
multi.matrix4 <- matrix[,1,drop=FALSE]
multi.matrix4 <- multi.matrix4[,j,drop=FALSE]
class(multi.matrix4)
}

# Deficiency of high trait values
for (j in 1 : N)
{matrix <- matrix(runif(total.species*2),total.species)
    class(matrix)
X <- 1 : total.species
X <- t(X)
# matrix[1:4,1] <- c(exp((X/22.6)^3)-1)
matrix[1:4,1] <- c(exp((X/13.55)^3)-1)
class(matrix)
matrix <- matrix[order(matrix( ,2))]
class(matrix)
# multi.matrix4[,j,drop=FALSE] <- matrix[,1,drop=FALSE]
# multi.matrix4[,j,drop=FALSE] <- matrix[,1]
# multi.matrix4[,j] <- matrix[,1,drop=FALSE]
# class(multi.matrix4)
multi.matrix4 <- matrix[,1,drop=FALSE]
multi.matrix4 <- multi.matrix4[,j,drop=FALSE]
class(multi.matrix4)
}

# Reading empirical data
Exp_07_2 <- file(description = "Exp_07_2", open = "r", blocking = TRUE, encoding = getOption("encoding"), raw = FALSE)
Exp_07_2 <- matrix(scan(Exp_07_2),nrow=4,byrow=TRUE)
read.matrix <- function(Exp_07_2){
    as.matrix(read.table(Exp_07_2))
}
Exp_07_2
class(Exp_07_2)
multi.matrix <- matrix(c(Exp_07_2),ncol=3)
class(multi.matrix)
multi.matrix <- multi.matrix(1:total.species,1:N)  
class(multi.matrix)
random <- runif(N)
multi.matrix2 <- t(multi.matrix)[order(t(multi.matrix)[,1], t(multi.matrix)[,2], t(multi.matrix)[,3], t(multi.matrix)[,4]),]
class(multi.matrix2) 
multi.matrix3 <- t(multi.matrix2)
class(multi.matrix3)
multi.matrix4 <- multi.matrix3[1:total.species,,drop=FALSE]
class(multi.matrix4)


# Adding a sorting column
random2 <- runif(total.species,1)
random2 <- multi.matrix4[,N+1,drop=FALSE]
sort.multi.matrix <- multi.matrix4(order(multi.matrix4[,1], multi.matrix4[,2], multi.matrix4[,3],multi.matrix4[,4]),N+1,drop=FALSE)

# loop adding one function at a time
for (j in 1 : J)

{loss.matrix <- sort.multi.matrix[nrow=species,ncol=j,drop=FALSE]
    class(loss.matrix)
max.value <- loss.matrix >= t
c(B) <- any(t(max.value),2)
c(C) <- all(c(B))
result(i,j) <- c(sum(C))
}
}

# Reporting
res <- mean(result)
res
t(res)

- user1842171

---

你可以选择将其翻译为“仅限于R”，或者使用Armadillo进行翻译，通过RcppArmadillo在R中使用。Armadillo的设计目标之一是使Matlab用户能够轻松进行此类转换。我已经成功地将其应用于一个昂贵的模拟问题，并获得了非常显著的速度提升。 - Dirk Eddelbuettel

也许如果您发布了带有“for”循环的R代码，我们可以更好地帮助您。只是猜测，但我怀疑multi.matrix4中只有N列，并且当j达到N+1时，循环会失败。 - nograpes

是的，没错。现在R代码已经发布。 - user1842171

Answer 1

虽然我手头没有Matlab和R，但我怀疑这就是问题的原因：

在R中，您试图为矩阵中不存在的位置分配值，结果：失败。

在Matlab中，您尝试为矩阵中不存在的位置分配值，结果：它会容忍您奇怪的选择、扩展您的矩阵并成功执行。

假设这就是问题，解决方法很简单：

在R中创建矩阵时，请确保它足够大，可以容纳后续要添加的所有内容。

这被称为初始化，在大多数情况下都是最佳实践。即使在Matlab中，也通常建议提前正确初始化变量，而不是让它们随着过程增长。

- Dennis Jaheruddin

---

谢谢。我将初始化代码，以便矩阵处于正确的比例。这很有意义。 - user1842171