2012-12-03 84 views
3

matlab中的代碼是爲了使生態系統中的物種在生態系統中丟失而發揮作用的可能性而創建的。現在,這段代碼必須翻譯成R.但是我有翻譯matlab中進行的矩陣操作的問題。將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) 
+0

你可以翻譯成'just R',或者你可以翻譯成[Armadillo](http://arma.sf.net),它可以通過[RcppArmadillo](http://dirk.eddelbuettel)從R中使用。 COM /代碼/ rcpp.armadillo.html)。 Armadillo的設計目標之一是爲Matlab用戶提供這種類型的轉換。我已經完成了昂貴的仿真問題,取得了巨大的成功,並且獲得了非常顯着的速度提升。 –

+0

也許如果你用'for'循環發佈R代碼,我們可以幫助你更好一點。只是猜測,但我懷疑'multi.matrix4'中只有'N'列,並且當'j'命中'N + 1'時循環失敗。 – nograpes

+0

是的,這是正確的。現在,R代碼已發佈。 – user1842171

回答

0

雖然我沒有Matlab的和R在手,讓我懷疑這是什麼造成的問題:

在R您嘗試分配給矩陣中不存在的位置,結果:它失敗

在Matlab中,您嘗試將其分配給矩陣中不存在的位置,結果是:它會原諒您的奇怪選擇,擴展您的矩陣併成功。

假設這就是問題所在,解決方法很簡單:

當創建R中的矩陣,確保它是足夠大的 包含所有你想在以後添加到它的東西。

這被稱爲initalization,並且在大多數情況下是最佳實踐。即使在Matlab中,通常建議儘可能提前初始化變量,然後讓它們隨着時間增長。

+0

謝謝。然後我將初始化代碼,使矩陣成正比。這說得通。 – user1842171