2013-02-21 89 views
0

我没有太多的经验,在R.我想写一个Gibbs采样在那里我有一个循环是这样的:如何向量化这个循环中的R

for (iNum in 1:totNum) { 
    rateNum <- Y3[iNum] 
    if(Y3[iNum] > 0) { 
     yStar3[iNum] <- rtnorm(1, mean = Mean3[iNum], sd = sqrt(Var3), 
       lower = gz[rateNum], upper = gz[rateNum + 1]) 
    } else if(Y3[iNum] == 0) { 
    yStar3[iNum] <- rtnorm(1, mean = Mean3[iNum], sd = sqrt(Var3), 
        lower = -Inf, upper = Inf); 
    } 
} 

这是服用过多的时间。我试图使用lapply,但这还不够快。有没有一种方法来实现这个循环的矢量化?

谢谢你,最好的问候!

+1

我会说是的,因为'rtnorm'被矢量,但你似乎对每个值使用不同的截断值(上,下),而且我认为这个参数没有被矢量化。 – joran 2013-02-21 19:54:26

回答

0

这是一个没有你的变量值的问题,但你想要做的事情是非常简单的。你想坚持所有的矢量化语句,并尽量不要占用过多的内存。这是基本策略:

第1步:找出如何计算所有数字。

# The number of values you need from 'rtnorm' 
sum(Y3 > 0) 
sum(Y3 == 0) 

# The means you need from the 'Mean3' array 
Mean3[Y3 > 0] 
Mean3[Y3 == 0] 

# Lower and upper limits for Y3 > 0 
gz[Y3[Y3 > 0]] 
gz[Y3[Y3 > 0] + 1] 

步骤2:上yStar3的向量滤波器使用这些值。我不能绝对肯定我所有的语法是完美的没有一些样本数据和变量值,但它应该是这个样子:

yStar3[Y3 > 0] <- rtnorm(
    sum(Y3 > 0), 
    mean = Mean3[Y3 > 0], 
    sd = sqrt(Var3), 
    lower = gz[Y3[Y3 > 0]], 
    upper = gz[Y3[Y3 > 0] + 1]) 

yStar3[Y3 == 0] <- rtnorm(
    sum(Y3 == 0), 
    mean = Mean3[Y3 == 0], 
    sd = sqrt(Var3), 
    lower = -Inf, 
    upper = Inf) 
+0

感谢大家的好评!我已经能够得到它的工作,并显着提高速度。这一步的一次迭代现在需要0.05秒而不是1秒。5秒,实际上为整个算法节省了数小时的计算时间! – user2096864 2013-02-22 10:18:52

+0

很高兴我们可以提供帮助。正如你可以说的那样,优化问题对我们很多人来说实际上是非常有趣的。下一次,如果您可以提供一些示例数据并要求优化而不是向量化,那么您可能会对您的问题产生更多兴趣和活动。 – Dinre 2013-02-22 18:19:23

+0

此外,如果您已达到解决方案,请继续选择您最满意的答案,或者通过单击其中一个答案旁边的复选标记回答您的原始问题最有帮助。这将有助于未来找到此页面的其他人并将结束该问题。 – Dinre 2013-02-25 12:36:36

2

因此,它不会出现,你必须迭代间的依赖性,这使得它非常简单的向量化

lhs = rtnorm(length(Y3), mean = Mean3, sd = sqrt(Var3), lower = gz[Y3], 
       upper = gz[Y3 + 1]) 
    rhs = rtnorm(length(Y3), mean=Mean3, sd = sqrt(Var3), lower=-Inf, upper=Inf) 

    ifelse(Y3 > 0, lhs, rhs) 

这里的问题是,rtnorm已经超过它的输入参数进行矢量化,意味着,下,和上。情况可能并非如此,在这种情况下,你将不得不做更多的工作。

+0

打败我吧。我相信'rtnorm'是msm包中的一个截断正态分布,所以它应该已经被矢量化了。 – user295691 2013-02-21 20:04:10

2

最简单的方法是生成两个条件的一半,并选择你想要的。该mean参数将采取载体手段,使你得到的东西是这样的:

yStar3 <- ifelse(
    Y3 > 0, 
    rtnorm(totNum, mean=Mean3, sd=sqrt(Var3), lower=gz[ratenum], upper=gz[rateNum+1]), 
    rtnorm(totNum, mean=Mean3, sd=sqrt(Var3), lower=-Inf, upper=Inf)) 

你必须与其中Y3小于零的情况下附加条件细化ifelse,潜在的,但是这是总体思路。

更新:@hadley表明移动ifelse的rtnorm内:

yStar3 <- rtnorm(totNum, mean=Mean3, sd=sqrt(Var3), 
    lower=ifelse(Y3>0,gz[rateNum], -Inf), 
    upper=ifelse(Y3>0,gz[rateNum+1], Inf)) 

现在有基本上是零不必要的计算回事。

更新:当然,1是错误的,正如评论者指出的;它应该是totNum。

+1

'ifelse'必须完全计算两个语句;你最好在'rtnorm'里面使用它' – hadley 2013-02-21 22:39:47

+0

你只是想到了我。 – user295691 2013-02-21 22:43:40

+1

您是否需要为'rtnorm'功能指定除'1'以外的数字?第一个参数告诉它返回多少个值。 – Dinre 2013-02-22 01:19:23