Skip to content

Instantly share code, notes, and snippets.

@hadley
Created Nov 19, 2012
Embed
What would you like to do?
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
double vacc3a(double age, bool female, double ily){
double p = 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily;
p = p * (female ? 1.25 : 0.75);
p = std::max(p, 0.0);
p = std::min(p, 1.0);
return p;
}
// [[Rcpp::export]]
NumericVector vacc3(NumericVector age, LogicalVector female, NumericVector ily) {
int n = age.size();
NumericVector out(n);
for(int i = 0; i < n; ++i) {
out[i] = vacc3a(age[i], female[i], ily[i]);
}
return out;
}
// [[Rcpp::export]]
NumericVector vacc4(NumericVector age, LogicalVector female, NumericVector ily) {
NumericVector p(age.size());
p = 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily;
p = p * ifelse(female, 1.25, 0.75);
p = pmax(0,p);
p = pmin(1,p);
return p;
}
// [[Rcpp::export]]
NumericVector vacc5(NumericVector age, LogicalVector female, NumericVector ily) {
NumericVector p(age.size());
p = 0.25 + 0.3 * 1 / (1 - exp(0.04 * noNA(age))) + 0.1 * noNA(ily);
p = noNA(p) * ifelse(noNA(female), 1.25, 0.75);
p = pmax(0, noNA(p));
p = pmin(1, noNA(p));
return p;
}
library(Rcpp)
library(microbenchmark)
vacc1a <- function(age, female, ily) {
p <- 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily
p <- p * if (female) 1.25 else 0.75
p <- max(0, p)
p <- min(1, p)
p
}
vacc1 <- function(age, female, ily) {
n <- length(age)
out <- numeric(n)
for (i in seq_len(n)) {
out[i] <- vacc1a(age[i], female[i], ily[i])
}
out
}
vacc2 <- function(age, female, ily) {
p <- 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily
p <- p * ifelse(female, 1.25, 0.75)
p <- pmax(0, p)
p <- pmin(1, p)
p
}
sourceCpp("vaccinate.cpp")
n <- 1000
age <- rnorm(n, mean = 50, sd = 10)
female <- sample(c(T, F), n, rep = TRUE)
ily <- sample(c(T, F), n, prob = c(0.8, 0.2), rep = TRUE)
stopifnot(
all.equal(vacc1(age, female, ily), vacc2(age, female, ily)),
all.equal(vacc1(age, female, ily), vacc3(age, female, ily)),
all.equal(vacc1(age, female, ily), vacc4(age, female, ily))
)
print(microbenchmark(
vacc1(age, female, ily),
vacc2(age, female, ily),
vacc3(age, female, ily),
vacc4(age, female, ily),
vacc5(age, female, ily)
))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment