Density

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector dpareto_cpp(const NumericVector& vec, const NumericVector& alpha, int x_m) {
  NumericVector ret(vec.length());
  for(auto i=0; i < vec.length(); ++i) {
    if (alpha[i] <= 0 || x_m <= 0) {
        ret[i] = NA_REAL;
    } else {
        ret[i] = alpha[i] * pow(x_m, alpha[i]) / pow(vec[i], alpha[i] + 1);
    }
  }   
  return ret;
}

dpareto_unknown <- function(vector, alpha, x_m) {
    density <- alpha * (x_m^alpha) / vector^(alpha + 1)
    density[vector < x_m] <- 0
    density[alpha <= 0] <- NaN
    density[x_m <= 0] <- NaN
    
    if (any(is.na(density))) {
        warning("All parameters should be positive.\n NaN is generated.\n")
    }

    return(density)
}

dpareto <- function(x, alpha, x_m){
    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)] <- 0
    
    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)] <- alpha_s * (x_m_s ^ alpha_s) / x_s^(alpha_s + 1)
    
    density
}

library(microbenchmark)
microbenchmark(
    dpareto_unknown(1:10000, 1:10000, 3),
    dpareto(1:10000, 1:10000, 3),
    dpareto_cpp(1:10000, 1:10000, 3)
)

## Unit: microseconds
##                                  expr       min        lq       mean
##  dpareto_unknown(1:10000, 1:10000, 3)   825.177   859.850  1055.2710
##          dpareto(1:10000, 1:10000, 3) 18225.937 18906.285 20233.4117
##      dpareto_cpp(1:10000, 1:10000, 3)   626.680   635.448   715.9069
##      median         uq       max neval
##    903.8795   928.3220  4470.585   100
##  19622.7630 20120.0480 48363.255   100
##    647.5475   661.2625  2064.273   100