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){

  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)) %>%
     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.


1 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)

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?

  • \$\begingroup\$ Thanks, all very good points! Will add profvis to the toolbox, looks promising. For larger dimensions, Stirling2 seems to be the largest bottleneck, as can be expected. About the mpfr - will have to test when can I really drop it; it was added to deal with extreme numeric values, both explosive and very close to 0, but maybe some approximations at both ends will work too. \$\endgroup\$
    – runr
    Apr 7, 2020 at 11:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.