Function for calculating Archimedean density functions
Some suggestions:


split code more, line by line, so it is easier to profile and see bottlenecks
use profvis
use %&gt;% 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 &lt;- function(x, tau, d){
  alpha &lt;- 1/tau
  l &lt;- sapply(1:d, function(k){
    p1 &lt;- Stirling2(n = d, k = k)
    p1 &lt;- as.numeric(p1)
    p1 * gamma(k - alpha) / gamma(1 - alpha) * x^(k-1)
  })
  sum(l)
}

joe_density &lt;- function(x, tau){
  d &lt;- length(x)
  h &lt;- prod((1 - (1 - x)^tau))

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

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


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