Goal

Write an R function to generate probabilities modeled by the following equation (IRT; 2 Parameter Logistic Model):

Data

set.seed(1)
a = runif(10, 1.5, 3)
b = rnorm(10, 0, 1)
theta = rnorm(10000000)
# the output of the implementation should result in
# a matrix with 100 mil. rows and 10 columns


Implementation

Strategy A

# the implementation of the equation
computeProbability = function(theta, b, a) {
x = exp(1) ^ (a * (theta - b))
return(x / (1 + x))
}

strategy.A = function(theta, b, a) {
n.rows = length(theta)
n.cols = length(b)

prob_mtx = matrix(nrow = n.rows, ncol = n.cols)

for(i in 1:n.rows)
{
prob_mtx[i, ] = computeProbability(theta[i], b, a)
}
return(prob_mtx)
}


Strategy B

strategy.B = function(theta, b, a) {
return(t(sapply(theta, computeProbability, b = b, a = a)))
}


Strategy C

strategy.C = function(theta, b, a) {
return(1 / (1 + exp(-sweep(outer(theta, b, "-"), 2, a, "*"))))
}


Timings

    # Strategy A        |       # Strategy B        |       # Strategy C
|                           |
user  system elapsed   |    user  system elapsed   |   user  system elapsed
64.76    0.27   65.08   |   82.01    0.91   82.93   |   7.81    0.64    8.46


Question: Strategy C is by far the most efficient way, but how can I make it even faster?

• That specific question isn't on-topic on this site, but we can review the code as it is instead. – Jamal Feb 21 '17 at 20:27
• Oh, my apologies. I thought it was too vague for StackOverflow. – Mihai Feb 21 '17 at 20:28
• Here is a great article that explains what is meant by vectorization in R. – jsb Feb 27 '17 at 16:31
• @Samuel It is really helpful. Thank you very much! – Mihai Feb 28 '17 at 15:36

Your strategy C is probably already nearly optimal. Here is an Rcpp based solution to squeeze out some additional seconds

library('Rcpp')
cppFunction(includes = '#include <math.h>',
code = 'Rcpp::NumericMatrix strategy_rcpp(const Rcpp::NumericVector& theta,
const Rcpp::NumericVector& b,
const Rcpp::NumericVector& a) {
const int n = theta.size();
const int m = a.size();
if (m != b.size())
Rcpp::stop("a and b must have equal length");
const double e = exp(1.0);
Rcpp::NumericMatrix ret(n, m);
for (int j = 0; j < m; ++j) {
for (int i = 0; i < n; ++i) {
ret(i, j) = 1.0 / (1.0 + pow(e, -a[j] * (theta[i] - b[j])));
}
}
return ret;
}')


library('microbenchmark')
microbenchmark(times = 1,
strategy.C(theta, b, a),
strategy_rcpp(theta, b, a)
)

Unit: seconds
expr       min        lq      mean    median        uq
strategy.C(theta, b, a) 10.996822 10.996822 10.996822 10.996822 10.996822
strategy_rcpp(theta, b, a)  6.921705  6.921705  6.921705  6.921705  6.921705