北京大学生物信息平台论坛

 找回密码
 立即注册
搜索
热搜: 通知 活动

【R语言】用向量和矩阵运算代替for循环,加快运算速度几...

[复制链接]
licheng 发表于 2016-5-20 00:07:08 | 显示全部楼层 |阅读模式
R语言中的for循环速度比较慢,apply系列函数也是循环。一个替代的办法是使用向量和矩阵运算,它们是R语言相对C语言的特色,虽然转换为二进制代码后也是循环,但它们在R语言中使用了特殊处理,运行的速度很快。

如果对全基因组的每个基因做t test,使用矩阵和向量自己计算所有基因的t统计量和p-value,要比调用t.test函数2万次快几十倍。
IMG_20160519_235543.png
IMG_20160519_235559.png
先看用牛顿方法递归寻找多项式根的例子:
IMG_20160519_235454.png
对根附近某一个起始值循环直到收敛,我可以得到一个根的估计:

  1. ############ Example: Newton's method for root finding

  2. x0 <- 1 #starting value
  3. x <- x0
  4. f <- x^3 + 2 * x^2 - 7
  5. tolerance <- 0.000001

  6. while (abs(f) > tolerance) {
  7.   f.prime <- 3 * x^2 + 4 * x
  8.   x <- x - f / f.prime
  9.   f <- x^3 + 2 * x^2 - 7
  10.   print(x)
  11. }

  12. x

  13. # confirm with plot
  14. curve(x^3 + 2 * x^2 - 7, -5, 5)
  15. abline(h = 0, col = "blue")
  16. abline(v = x, col = "red")
复制代码
这段代码可以几乎原封不动地改为对向量x做计算,因为R语言可以自动识别x是一个数还是一个向量,并做相应运算。这样我们就可以对一组起始值同时循环牛顿方法,检验它们是否都能收敛到对根的估计,这样就避免了双重循环:

  1. ### starting values as a vector
  2. x0 <- seq(-3, 3, length = 30)

  3. x <- x0
  4. f <- x^3 + 2 * x^2 - 7
  5. tolerance <- 0.000001

  6. while (max(abs(f)) > tolerance) {   # note the change
  7.   f.prime <- 3 * x^2 + 4 * x
  8.   x <- x - f / f.prime
  9.   f <- x^3 + 2 * x^2 - 7
  10.   #print(x)
  11. }

  12. x
复制代码
现在到我们的挑战任务了,你能否先写出R代码,对两个向量计算出一个t检验统计量,然后扩展代码,对两个矩阵进行向量和矩阵操作,计算出对应所有基因(行)的t统计量的向量?
IMG_20160519_235333.png
参考这个链接的课件和R代码:
https://yunpan.cn/cSFFUesvRXmN7

回复

使用道具 举报

北京大学生物信息平台论坛

GMT+8, 2017-11-19 22:11 , Processed in 0.085974 second(s), 26 queries .

Powered by Discuz! X3

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表