Re: [問題] grid生樣本 避免for迴圈

[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/
: ※ 引述《ntpuisbest (阿龍)》之銘言：
: : n <- 20
: : theta=seq(0,1,len=100)
: : rt(n ,1/theta )
: : 如題
: : 我想要生100組 ，每組都是n=20的t分配樣本
: : 只是這100組的theta都不一樣
: : 我像上面那樣打 只會回傳20個樣本
: : 並不是我想要的 2000個樣本 請問要如何打才能要我要的結果
: : 想避免for loop
: : 用loop的話 我知道怎麼做

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 119.14.59.166
※ 文章網址: https://www.ptt.cc/bbs/R_Language/M.1553702206.A.0CA.html
※ 編輯: celestialgod (119.14.59.166), 03/28/2019 00:10:33
推 locka: 竟然…！(登愣XDD
→ locka: 但還是想請問c大：1. 當初怎麼知道Vectorize 背後用的是 ma
→ locka: pply？
乖乖在console裡面輸入mapply就可以看到了XD

→ locka: 2. 如何判斷 rt() 沒辦法用向量化的方式做？
→ locka: 還是必須說，大大太神了
看文件，df不支援放向量

推 locka: 自問自答1：help("Vectorize")就有了。
→ locka: 衍生問題3：所以 Vectorize(f) 只是單純把又臭又長的mapply
→ locka: (f…)包起來，實際上跟直接用mapply比並沒有差？
Vectorize只是幫你用mapply包起來沒錯，比直接用mapply還慢

因為會有一些overhead...

建議就是直接用迴圈包起來，像我上面那樣做就好XD

推 locka: df degrees of freedom (> 0, maybe non-integer). df = In
→ locka: f is allowed.
→ locka: 不是只是說 df 參數可以是無限大(Inf)的意思嗎？還是我理解
對 (Inf 就normal)

→ locka: 錯誤？
→ locka: 哪裡有說他不能用vector of degrees of freedom？
基本上沒寫可以就是不行XDD，不然像是n, prob什麼都會寫可以輸入vector
※ 編輯: celestialgod (119.14.59.166), 03/28/2019 00:39:10
推 locka: 感謝釋疑，但又有疑問了(舉手) 最後一個 >"<
→ locka: 所以一般所謂R的向量化到底是什麼意思啊？？如果*apply, Ve
→ locka: ctorize 都不是的話。(先謝謝版主大大一一耐心回覆！)
我覺得這篇寫得很好，先看這篇吧XD

http://alyssafrazee.com/2014/01/29/vectorization.html
※ 編輯: celestialgod (118.163.170.73), 03/28/2019 14:00:37
→ andrew43: http://bit.ly/2WqRtYt 這篇大概把幾種（但不限於R）向
→ andrew43: 量化的特徵、速度快的（多種）原因、以及何時還是用for
→ andrew43: 就好。一點點資工知識應該就可以看懂。
→ locka: 樓上這篇跟我上一篇貼的一樣啊XD
推 andrew43: 喔。XD。
→ obarisk: 新版本的R應該不太需要用compiler了。
→ celestialgod: 基於測試，還是加進去看看