Function for calculating Archimedean density functions

Question

I am not sure whether this is the right place to ask this, so please correct me if it's not.

I am following this paper (e.g., Corollary 3.3) for some derived functional representations of the densities and writing code for those. The code for one of the families:

(x will be a vector of 1xd between [0,1], tau - scalar in [1, infty).)

'joe_density' <- function(x, tau){
  require(gmp)
  require(Rmpfr)
  require(magrittr)

  d <- length(x)

  P_d <- function(x, tau, d){
    alpha <- 1/tau
    lapply(1:d, function(k){
      as.numeric(gmp::Stirling2(n = d, k = k)) * gamma(k - alpha) / gamma(1 - alpha) * x^(k-1)
    }) %>% do.call(sum, .)
  }

  h <- prod((1 - (1 - x)^tau))

  mpfr( tau^(d-1), 16) * 
    mpfr( prod((1 - x)^(tau - 1)), 16) / 
    mpfr( (1 - h)^(1 - 1/tau), 16) * 
    mpfr( P_d( h/(1-h), tau, d), 16)
}

Example call:

# Simple case
x <- pnorm(rnorm(100))
joe_density(x, tau = 2)

# Longer performance time when looping, d=2
expand.grid(x = seq(-5,5,by=0.1), y = seq(-5,5,by=0.1)) %>%
  apply(.,1,function(z){
     x <- pnorm(z)
     joe_density(x, tau = 2)
  })

In essence, the typical calculations require tons of looping, gamma/Stirling numbers, and I have noticed that the default precision is usually not enough, hence the Rmpfr usage. One time evaluation of the code is not a problem, however, the calculations become very lengthy when the function is called thousands of time (e.g., passing it to optim and the likes).

Are there simpler ways to program this? This particular function is just an example, but most of the different functions follow similar approach, with a high usage of Stirling/gamma, lengthy products and etc.

Could, for example, rewriting parts of the code in Rcpp be helpful, or is this not the type of task that we could expect gains here? I have no experience in writing cpp code, so am not sure what to expect.

Answer 1

Some suggestions:

• split code more, line by line, so it is easier to profile and see bottlenecks

• use profvis

• use %>% less often if that part of code is called frequently, it has some overhead

• separate function definitions

This should run a little bit faster: (15s vs 27s on your example)

P_d <- function(x, tau, d){
  alpha <- 1/tau
  l <- sapply(1:d, function(k){
    p1 <- Stirling2(n = d, k = k)
    p1 <- as.numeric(p1)
    p1 * gamma(k - alpha) / gamma(1 - alpha) * x^(k-1)
  })
  sum(l)
}

joe_density <- function(x, tau){
  d <- length(x)
  h <- prod((1 - (1 - x)^tau))

  v1 <- tau^(d-1)
  v2 <- prod((1 - x)^(tau - 1))
  v3 <- (1 - h)^(1 - 1/tau)
  v4 <- P_d( h/(1-h), tau, d)

  p1 <- mpfr(v1, 16) # slow
  p1 * v2 / v3 * v4 # slow
}

The slowest part is mpfr & the last line, because each value is converted to mpfr. Maybe you can deal with precision loss in some kind of different way?