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.