Re: [问题] grid生样本 避免for循环

楼主: ntpuisbest (阿龙)   2019-04-07 10:51:34
※ 引述《celestialgod (天)》之铭言:
: 1. Vectorize用的mapply,其实还是循环,而且表现更糟
: 2. rt没有办法直接执行原PO想要的那种方式,你可以自己用Rcpp刻一个
: 3. 花时间查资料,有时候不如看文件或是source code来的有效
: Vectorize直接看应该不难找到重点
: rt可以直接看文件
: 会看到df degrees of freedom (> 0, maybe non-integer). df = Inf is allowed.
: 这里很明显告诉你他不能用vector of degrees of freedom....
: 下面我比较一下几种做法:
: n <- 20
: theta <- seq(0, 1, len=100) + 0.001 # add small value to avoid warning
: for_loop_rt <- function(n, theta) {
: stopifnot(length(n) == 1, abs(n - floor(n)) < 1e-8)
: X <- matrix(NA, n, length(theta))
: for (i in seq_along(theta))
: X[ , i] <- rt(n, theta[i])
: return(X)
: }
: library(compiler)
: for_loop_rt_compiled <- cmpfun(for_loop_rt)
: rt_vectorized <- Vectorize(rt)
: library(microbenchmark)
: microbenchmark(
: sapply = sapply(theta, function(t) rt(20, 1/t)),
: vectorized = rt_vectorized(rep(n, length(theta)), theta, ncp = rep(0,
: length(theta))),
: for_loop = for_loop_rt(n, theta),
: for_loop_cmpfun = for_loop_rt_compiled(n, theta)
: )
: # Unit: microseconds
: # expr min lq mean median uq max neval
: # sapply 556.8 576.10 597.278 584.30 594.25 1639.1 100
: # vectorized 898.6 923.45 1041.179 936.45 961.45 2637.4 100
: # for_loop 546.4 557.00 649.992 564.85 571.65 5844.1 100
: # for_loop_cmpfun 541.9 558.65 679.358 566.45 573.85 2290.1 100
: 这结果,很讽刺地告诉你用Vectorize只是让事情更糟而已XDDDDDD
: 然后反而用for loop最快.... 这件事情其实我有在板上谈过....
: matrix操作 R真的很快.... 试着相信R一下^.<
: ※ 引述《locka (locka)》之铭言:
: : 感谢 celestialgod 版主大大提点:
: : 以前以为 *apply 家族的函数就已经是向量化(vectorized)的写法了
: : 查了资料才发现其实底层背后还是有 for 循环 (觉得震撼啊...)
: : 试试看这样的写法
: : theta <- seq(0,1,len=100)
: : df <- rep(19,len=100)
: : n <- rep(20,len=100)
: : vrt <- Vectorize(rt)
: : x <- vrt(n=n, df=df, ncp=1/theta)
: : 于是 x[,1] ... x[,100] 就是100个 n 等于20 然后对应各自 delta 值的 t 分配样本了
: : (但是不知道 df, n 的预先定义有没有意义?)
: : 请版上各位高手再指点~ 谢谢大家
: : ======
: : 补充:
: : 但还是有查到 *apply function 的好处:
: : 1. 程式易读性
: : 2. 会 pre-allocate 向量的内存空间
: : 2. 只影响区域变量不会改变全域变量
: : ref: https://www.r-bloggers.com/vectorization-in-r-why/
感谢大家回答
其实最后我选择用第一篇贴文的lapply来做
我的叙述不太好
我并不是不想用循环
只是不想写for 循环
对于我很晚回复这件事我也感到抱歉
毕竟大家都好热烈的讨论
果然R版真的是藏龙卧虎XDD

Links booklink

Contact Us: admin [ a t ] ucptt.com