본문으로 바로가기



Make your own R package! - 6. R 함수 벡터화(Vectorized function)에 대하여

category R Programming/Rpackage 2017.08.02 02:05
Make your own R package! - 6. 함수 벡터화(Vectorized function)에 대하여

오늘의 주제는 함수를 정의하는 방법에 대한 내용이다. 패키지에 함수를 정의해서 사용하는 방법은 사실 이전 포스팅에서 많이 다루었다. 하지만 만약 훗날 다른 사람들을 위하여 함수를 정의하고 배포할 경우 고려해야하는 사항들을 한번 알아보도록 하자. 이 포스팅은 필자가 수강한 학교 수업의 숙제에서 따온 내용이다.

Pareto분포의 확률밀도함수 dpareto 정의하기

R은 여러 통계 분포에 관한 함수를 제공하고 있다. 예를 들어 정규분포의 경우, 확률밀도함수는 dnorm, 누적분포함수는 pnorm 등 앞의 분포의 이름 앞에 제공하는 함수의 성질을 간략한 약자를 사용하여 나타내고있다. 오늘 우리가 만들어볼 함수는 dpareto 함수, 즉 파레토 분포의 확률밀도함수이다. 위키피디아의 파레토 분포의 PDF 정의를 살펴보면 다음과 같다.

알파에 따른 파레토 확률밀도함수 형태(사진출처: 위키피디아)

\[ \frac{\alpha\,x_\mathrm{m}^\alpha}{x^{\alpha+1}}\text{ for }x\ge x_m \]

여기서 \(\alpha > 0\)는 실수값을 갖는 shape 모수, \(x_m > 0\)의 경우는 실수값을 갖는 scale 모수이다. 함수의 정의역은 \([x_m, +\infty)\)이며, 정의역 밖의 함수값은 밀도함수의 경우 0으로 처리한다.

위의 함수를 패키지에 넣어서 배포를 한다고 가정하자. R을 처음 사용하는 독자의 경우는 다음과 같이 정의할 수 있다.

paretodens <- function(x, alpha, x_m){
  if (x < x_m) {
    density <- 0
  } else {
    density <- alpha * (x_m ^ alpha) / x^(alpha + 1)  
  }
  density
}

위의 paretodens 함수는 혼자서 사용하기에는 좋다. 왜냐하면 사용자가 자기 자신이고, 이미 이 함수를 짠 시점에서 파레토 확률밀도함수에 대한 사전지식을 가지고 있기 때문에, shape과 scale이 어떠한 값을 가져야하는지 알 것이고, 정의역이 어떻게 정의되었는지 알고 있을 것이기 때문이다. 따라서 위의 정의된 함수를 잘못 사용할 가능성이 적다. 하지만 이것을 다른 독자들에게 배포를 한다고 생각하면 생각해야할 것들이 몇가지가 존재한다.

잘못 입력한 모수에 대한 에러 출력

사용자들 중에는 알맞는 모수가 어떤지에대한 사전 지식이 없는 사용자가 있을 것이다. 이들을 위하여 잘못 입력을 한 모수들이 들어왔을 경우 다음과 같이 에러 처리를 해주는 것이 좋다.

paretodens <- function(x, alpha, x_m){
  if (alpha < 0 || x_m < 0){
    warning("All parameters should be positive.")
    return(NaN) 
  }
  if (x < x_m) {
    density <- 0
  } else {
    density <- alpha * (x_m ^ alpha) / x^(alpha + 1)  
  }
  density
}

density를 NaN으로 처리해주는 것은 다른 정의된 함수들의 에러처리를 따르도록 하기 위해서이다. R에서 정의된 통계분포 함수들은 모수의 조건을 위반했을 경우 NaN으로 에러처리를 하고있다.

dgamma(3, -1, 1)
## Warning in dgamma(3, -1, 1): NaNs produced
## [1] NaN

log 기능 고려하기

현재 대부분의 통계 함수들은 log 옵션을 가지고 있다. 로그 옵션을 부여하는 것은 사용자의 편의도 고려할 수 있을 뿐 아니라 계산 과정에서의 소수점 부분의 손실도 줄일 수 있기 때문에 로그를 취한 후 계산을 하고 마지막에 exp 함수를 사용하여 원래대로 만들어주는 것이 좋다.

paretodens1 <- function(x, alpha, x_m, log = FALSE){
  if (alpha < 0 || x_m < 0){
    warning("All parameters should be positive.")
    return(NaN) 
  }

  if (x < x_m) {
    density <- -Inf
  } else {
    density <- log(alpha) + alpha * log(x_m) - (alpha + 1) * log(x)
  }

  if (log) density
  else exp(density)
}
paretodens(5, 2, 3)
## [1] 0.144
paretodens1(5, 2, 3)
## [1] 0.144
paretodens1(5, 2, 3, log = TRUE)
## [1] -1.937942

벡터연산 가능 함수 만들기

R 프로그래밍의 묘미는 코딩을 얼마나 벡터화 시키는가에 있다고해도 과언이 아니다. 따라서 우리가 배포할 함수들을 사용하는 사용자들은 당연히(?) 혹은 개념있는(?) 패키지라면 함수들이 벡터화된 입력값을 잘 처리할 수 있을 것이라 기대한다. 따라서 다음과 같은 명령어들을 처리할 수 있도록 해줘야 한다.

dpareto(1:3, 2, 3)
dpareto(3, 1:3, 3)

위의 입력값들을 처리하기 위해서 다음과 같이 입력값을 매트릭스 형식으로 바꿔준다.

x <- 4:1
alpha <- -2:5
x_m <- 1:3 

input <- paste(x, alpha, x_m)
input <- as.numeric(unlist(strsplit(input," ")))
input <- matrix(input, ncol=3, byrow=T)
input
##      [,1] [,2] [,3]
## [1,]    4   -2    1
## [2,]    3   -1    2
## [3,]    2    0    3
## [4,]    1    1    1
## [5,]    4    2    2
## [6,]    3    3    3
## [7,]    2    4    1
## [8,]    1    5    2

위의 테크닉을 사용해서 dpareto 함수를 정의하고, 앞에서 정의한 paretodens함수와 연동시키자.

dpareto1 <- function(x, alpha, x_m, log = FALSE){
  input <- paste(x, alpha, x_m)
  input <- as.numeric(unlist(strsplit(input," ")))
  input <- matrix(input, ncol=3, byrow=T)
  density <- apply(input, 1, paretodens2)
  if (any(is.nan(density))) {
        warning("All parameters should be positive.\n NaN is generated.\n")
  }
  if (log) density
  else exp(density)
}
paretodens2 <- function(input){
  x <- input[1]
  alpha <- input[2]
  x_m <- input[3]

  if (alpha <= 0 | x_m <= 0){
    value <- NaN
  } else if (x < x_m) {
    value <- -Inf
  } else {
    value <- log(alpha) + alpha * log(x_m) - (alpha + 1) * log(x)
  }
  value
}

위의 dpareto1 함수는 다음과 같은 인풋을 처리할 수 있게된다.

dpareto1(1:4, 2, 3)
## [1] 0.0000000 0.0000000 0.6666667 0.2812500
dpareto1(3, -2:4, 3)
## Warning in dpareto1(3, -2:4, 3): All parameters should be positive.
##  NaN is generated.
## [1]       NaN       NaN       NaN 0.3333333 0.6666667 1.0000000 1.3333333
dpareto1(4, -2:4, 1:7)
## Warning in dpareto1(4, -2:4, 1:7): All parameters should be positive.
##  NaN is generated.
## [1]  NaN  NaN  NaN 0.25 0.00 0.00 0.00

위의 두 함수를 패키지에 넣는다면 dpareto의 경우는 export를 가능하게 처리하고, paretodens의 경우는 사용자가 사용하지 못하도록 export를 금지시켜 놓는 것이 좋을 것이다.

순수 벡터 연산으로 속도 빠르게하기

위의 dpareto1의 경우, apply를 사용함으로서 코드가 간결해졌지만, 완전 벡터화를 시킨 코드보다는 시간이 많이 소요된다. 따라서 패키지에 들어가는 함수의 경우 다음과 같이 벡터화된 dpareto함수가 더 선호될 것이다.

Version 1

dpareto2 <- function(x, alpha, x_m, log = FALSE){
  input <- paste(x, alpha, x_m)
  input <- as.numeric(unlist(strsplit(input," ")))
  input <- matrix(input, ncol=3, byrow=T)
  
  x <- input[,1]
  alpha <- input[,2]
  x_m <- input[,3]
  
  density <- numeric(length(x))
  invalid <- alpha <= 0 | x_m <= 0
  
  density[invalid] <- NaN
  num <- sum(invalid)
  if (num > 0){
    warning("All parameters should be positive.\n NaN is generated.\n")
  }
  
  density[!invalid & (x < x_m)] <- -Inf
  
  alpha_s <- alpha[!invalid & !(x < x_m)]
  x_m_s <- x_m[!invalid & !(x < x_m)]
  x_s <- x[!invalid & !(x < x_m)]
  
  density[!invalid & !(x < x_m)] <- 
      log(alpha_s) + alpha_s * log(x_m_s) - (alpha_s + 1) * log(x_s)
  
  if (log) density
  else exp(density)
}

결과를 비교해보도록하자.

all.equal(dpareto1(1:4, 2, 3), dpareto2(1:4, 2, 3))
## [1] TRUE
all.equal(dpareto1(3, -2:4, 3), dpareto2(3, -2:4, 3))
## [1] TRUE
all.equal(dpareto1(4, -2:4, 1:7), dpareto2(4, -2:4, 1:7))
## [1] TRUE

완전 벡터형식으로 짠 dpareto2는 속도면에서 dpareto1보다 우월한 것을 확인할 수 있다.

library(microbenchmark)
set.seed(2017)
microbenchmark(
  dpareto1(1:10000, 1:10000, 3),
  dpareto2(1:10000, 1:10000, 3)
)
## Unit: milliseconds
##                           expr      min       lq     mean   median
##  dpareto1(1:10000, 1:10000, 3) 63.73826 67.61586 73.95496 72.00829
##  dpareto2(1:10000, 1:10000, 3) 25.63616 26.44239 30.35525 29.10025
##        uq       max neval
##  78.14867 138.05569   100
##  31.46446  88.92997   100

Version 2

위의 코드를 좀 더 다듬어 보면 다음과 같다. 앞의 input값을 만드는 구간이 병목 구간이었던 것을 바꿔서 보았다. 또한 version 2의 코드는 Version 1과 비교해서 읽기가 수월해졌다.

dpareto3 <- function(x, alpha, x_m, log = FALSE) {
  mlength <- max(length(x), length(alpha), length(x_m))

  density <- vector(mode="numeric", length = mlength )
  x     <- rep(x, length.out = mlength)
  alpha <- rep(alpha, length.out = mlength)
  x_m   <- rep(x_m, length.out = mlength)

  #Filterling error condition & support condition
  invalid <- alpha <= 0 | x_m <= 0
  sub <- x >= x_m & alpha > 0 & x_m > 0

  # Invalid parameter
  if (sum(invalid) > 0){
    warning("All parameters should be positive.\n NaN is generated.\n")
  }
  
  #Give 'NaN' values for the inputs which contains errors.
  density[sub] <- log(alpha[sub]) + alpha[sub] * log(x_m[sub]) -
                  (alpha[sub] + 1) * log(x[sub])
  density[invalid] <- NaN
  density[!invalid & !sub] <- -Inf

  if (log) density
  else exp(density)
}

dpareto2과 dpareto3 함수의 성능 비교해보면, 성능면에서도 dpareto2보다 dpareto3가 우월한 성능을 보이는 것을 알 수 있다.

all.equal(dpareto2(1:4, 2:4, 3), dpareto3(1:4, 2:4, 3))
## [1] TRUE
all.equal(dpareto2(3:10, -2:4, 3:5), dpareto3(3:10, -2:4, 3:5))
## [1] TRUE
all.equal(dpareto2(4, -2:4, 1:7), dpareto3(4, -2:4, 1:7))
## [1] TRUE
set.seed(2017)
microbenchmark(
  dpareto2(1:10000, 1:10000, 3),
  dpareto3(1:10000, 1:10000, 3)
)
## Unit: milliseconds
##                           expr       min        lq      mean    median
##  dpareto2(1:10000, 1:10000, 3) 25.692393 26.106560 28.619484 27.563698
##  dpareto3(1:10000, 1:10000, 3)  1.382407  1.402502  1.941973  1.426029
##         uq       max neval
##  27.996047 84.794685   100
##   3.041538  3.570056   100

Version 3

다음은 페이스북 R Korea - KRSG(Korean R Study Group) 그룹의 Kyung Mo Kweon님의 코딩이다.

dpareto4 <- function(x, alpha, x_m, log = FALSE) {
  density <- log(alpha) + alpha * log(x_m) - (alpha + 1) * log(x)
  density[x < x_m] <- -Inf
  density[alpha <= 0 | x_m <= 0] <- NaN
    
  if (any(is.nan(density))) {
    warning("All parameters should be positive.\n NaN is generated.\n")
  }

  if (log) density
  else exp(density)
}

앞부분의 input부분을 R에서 recycling을 이용하여 대체하고, 처음 density 벡터에 쓰레기값을 먼저 생성한 다음 코드 뒷부분에 정의역 밖에서의 0값과 \(\alpha\)\(x_m\)의 검정을 처리하는 것을 볼 수 있다.(이렇게 함으로서 필자가 처음에 짰던 조건문들을 다 피해서 짠 것을 볼 수 있다.) 테스트를 위하여 코드를 돌려보았다. 결과는 통과!

all.equal(dpareto3(1:4, 2:4, 3:8), dpareto4(1:4, 2:4, 3:8))
## [1] TRUE
all.equal(dpareto3(3, -2:4, 3:5), dpareto4(3, -2:4, 3:5))
## [1] TRUE
all.equal(dpareto3(4, -2:4, 1:7), dpareto4(4, -2:4, 1:7))
## [1] TRUE
stopifnot(all.equal(dpareto4(3,1,2), 0.2222222222))
stopifnot(all.equal(dpareto4(1,3,2), 0.0))
stopifnot(all.equal(dpareto4(1,-3,2), NaN))
stopifnot(all.equal(dpareto4(1, 3,-2), NaN))
stopifnot(all.equal(dpareto4(3:5,1, 2), c(0.2222222222, 0.1250000, 0.0800000)))
stopifnot(all.equal(dpareto4(6, 1, 2:4), c(0.05555555556, 0.08333333333, 0.11111111111)))
dpareto4(3, 1, 2, log = TRUE) # -1.504077
## [1] -1.504077

속도 역시 필자의 version2 보다 빠르고 max time부분에서 훨씬 안정적인 것을 확인했다. Rcpp로 쓴 함수와의 추가적인 비교는 다음의 링크를 참고하자.

set.seed(2017)
microbenchmark(
  dpareto1(1:1000, 1:10, 3:6),
  dpareto2(1:1000, 1:10, 3:6),
  dpareto3(1:1000, 1:10, 3:6),
  dpareto4(1:1000, 1:10, 3:6)
)
## Unit: microseconds
##                         expr      min        lq      mean   median
##  dpareto1(1:1000, 1:10, 3:6) 5524.000 5621.8165 6006.1433 5745.406
##  dpareto2(1:1000, 1:10, 3:6) 2165.163 2182.5970 2303.8600 2211.015
##  dpareto3(1:1000, 1:10, 3:6)  152.094  156.3765  165.2434  159.473
##  dpareto4(1:1000, 1:10, 3:6)  105.059  112.9055  135.0023  127.552
##        uq      max neval
##  5993.660 8953.825   100
##  2326.108 3713.143   100
##   164.406  294.782   100
##   150.227  255.395   100

Reference

[1] U of Iowa, Computer Intensive Statistics: http://homepage.divms.uiowa.edu/~luke/classes/STAT7400/


SHARE TO

신고


티스토리 툴바