2016-06-09 68 views
-1

我有一个简短的脚本,它用两个不同的分布进行采样来模拟物理过程。看评论。我怎样才能将迭代结果放入矩阵进行进一步的统计分析?我查看了以前回答的问题,但仍然无法使其工作。我明白,如果循环不是R中的首选方法,但循环是我基于基本的Perl和Python理解的,其他所有内容都会让我困惑。R新手:将一个循环的输出保存到一个矩阵中

library(truncnorm)          
library(mc2d)          
o <- 0.04          
n <- 10  # number of random samples - kept low for debugging 
md <- seq(0,0.70,by=0.05) # md for mode in the PERT distribution 
for(i in md) { # iterates over all modes in PERT distribution      
f <- rpert(n, min=0, mode=md, max=.99, shape=4) # samples from PERT distribution        
a <- rtruncnorm(n, a=0, b=Inf, mean = 5.44, sd = 0.43) # samples from normal distribution     

ma <- a*(1-f)+ f*o # calculates results       

print(ma) # I need this in a matrix 
} 
+1

忘记你从其他编程语言知道什么。相反,记住你的代数课并思考向量。如果你想为矢量的每一个元素做一些事情,很有可能有一个函数可以一次完成整个矢量。 – Roland

+0

嗨@Klaus,下面的答案是否回答你的问题?如果是这样,请随时通过点击旁边的“v”符号来接受您认为最有用的答案。另请参阅http://stackoverflow.com/tour如果您的问题尚未解决,请随时提供有关不起作用的说明。 – coffeinjunky

回答

0

你可以列所需数量的动态初始化一个空矩阵,然后rbind模拟结果到矩阵:

library(truncnorm)          
library(mc2d)          
o <- 0.04          
n <- 10  # number of random samples - kept low for debugging 
md <- seq(0,0.70,by=0.05) # md for mode in the PERT distribution 

resultMatrix <- matrix(nrow = 0, ncol = n) 
for(i in md) { # iterates over all modes in PERT distribution      
    f <- rpert(n, min=0, mode=md, max=.99, shape=4) # samples from PERT distribution        
    a <- rtruncnorm(n, a=0, b=Inf, mean = 5.44, sd = 0.43) # samples from normal distribution       
    ma <- a*(1-f)+ f*o # calculates results       
    resultMatrix <- rbind(resultMatrix, ma) 
} 

resultMatrix 
     [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8] 
ma 4.871186 4.527776 6.155239 2.687423 3.172824 4.150493 4.309858 4.058311 
ma 4.492796 4.023142 5.202100 3.549102 3.121166 4.030550 4.072121 4.754218 
ma 5.209793 3.604659 4.546211 4.019750 4.230491 2.230483 5.321552 4.194712 
ma 5.588713 5.845059 4.673996 5.068933 4.445595 2.556984 3.206710 2.198054 
ma 5.161478 4.224916 4.731809 4.337504 5.279807 4.800994 1.920244 2.287079 
ma 4.370172 4.174557 4.657961 3.551273 3.176608 1.967586 2.082636 2.382379 
ma 4.752645 5.025765 4.077100 4.035936 3.696961 4.376170 4.756229 3.799819 
ma 4.895087 4.476946 4.849364 4.948230 4.254575 2.860149 4.385260 4.909378 
ma 4.600699 2.955893 4.011127 5.018648 3.446420 2.684369 3.733717 3.784529 
ma 3.015752 5.150226 4.641076 5.403140 1.566149 4.467812 2.624535 4.140788 
ma 5.126098 4.311545 3.244769 2.922413 5.712901 2.981147 2.106302 4.173604 
ma 4.497303 5.128249 2.420177 3.460802 3.158532 1.826757 3.705091 2.092096 
ma 5.133561 3.937170 5.159742 3.097803 3.157485 4.583058 3.529645 5.299575 
ma 5.384795 5.977110 4.269142 3.898964 5.024477 2.174062 4.364693 2.060221 
ma 4.963641 4.670167 3.732576 5.075890 3.384682 3.581102 2.846963 4.388798 
     [,9] [,10] 
ma 2.035035 3.047999 
ma 3.119995 1.962640 
ma 5.107974 3.442291 
ma 2.810294 1.258326 
ma 2.951629 2.220695 
ma 4.461524 2.804621 
ma 3.077005 5.595100 
ma 1.248858 3.427024 
ma 1.983569 2.546225 
ma 1.510322 2.482890 
ma 4.652874 1.692841 
ma 4.029131 4.909566 
ma 3.502746 2.017282 
ma 1.992627 3.583658 
ma 3.336058 3.147302 
+3

这对于大型MD来说可能非常缓慢。动态增长的对象,例如'rbind()'通常是主要的性能瓶颈,通常不被认为是好的做法。我建议在循环之前在其定义处分配'resultMatrix'所需的维度,并使用循环内的行索引填充数据行。 – RHertel

0

你可以试试这个:

output <- matrix(

    unlist(

    lapply(md, function(md){ 

     f <- rpert(n, min=0, mode=md, max=.99, shape=4) # samples from PERT distribution        
     a <- rtruncnorm(n, a=0, b=Inf, mean = 5.44, sd = 0.43) # samples from normal distribution     
     ma <- a*(1-f)+ f*o # calculates results 

    }) 

), 

    ncol=10, byrow=TRUE 

) 
1

我无法看到循环的结果如何取决于i,因此可以避免直接使用矩阵而避免循环。

library(truncnorm); library(mc2d)          
o <- 0.04; n <- 10 
md <- seq(0,0.70,by=0.05) 

pertmat <- matrix(rpert(n*length(md), min=0, mode=md, max=.99, shape=4), 
        ncol=n, nrow=length(md)) 
amat <- matrix(rtruncnorm(n*length(md), a=0, b=Inf, mean = 5.44, sd = 0.43), 
       ncol=n, nrow=length(md))   

finalmat <- amat*(1-pertmat) + pertmat*o 

需要注意的是,此刻,你传递整个矢量mdrpert函数作为参数。如果您想要传递md的单个元素,则必须将md参数改为i。更明确地说,比较以下两个输出较短md矢量(为了避免混乱屏幕空间):

short_md <- seq(0,0.70,by=0.35) 
for(i in short_md) print(i) # single elements 
[1] 0 
[1] 0.35 
[1] 0.7 
for(i in short_md) print(short_md) # in each iteration, the entire vector is used 
[1] 0.00 0.35 0.70 
[1] 0.00 0.35 0.70 
[1] 0.00 0.35 0.70 

要想从这个矩阵,您可以使用:

pertmat_i <- t(sapply(md, function(x) rpert(n, min=0, mode=x, max=.99, shape=4))) 

,然后用这个上面的矩阵。