[REVIEW] Implementation of the Lévy-frailty model

[hsloot]

Checklist

• Rcpp implementation
• Original R implementation
• Bivariate R implementation

Proposal

Review the implementation of the exogenous shock model in Rcpp:

rmo/src/rmo.cpp

Lines 289 to 310 in cf85b82

	NumericMatrix Rcpp__rmo_lfm_cpp(unsigned int n, unsigned int d, double rate, double rate_killing, double rate_drift, Function rjump, List rjump_arg_list) {
	NumericVector unit_exponentials(d);
	NumericMatrix cpp_subordinator;
	unsigned int count;
	NumericMatrix out(n, d);
	for (unsigned int k=0; k<n; k++) {
	unit_exponentials = Rcpp::rexp(d);
	cpp_subordinator = sample_cpp(rate, rate_killing, rate_drift, rjump, rjump_arg_list, unit_exponentials);
	for (unsigned int i=0; i<d; i++) {
	count = 0;
	while (cpp_subordinator(count, 1) < unit_exponentials[i] && count < cpp_subordinator.nrow())
	count += 1;
	if (cpp_subordinator.nrow() == count)
	stop("internal error: exponential value out of subordinator range");
	out(k, i) = cpp_subordinator(count, 0);
	}
	}
	return out;
	}

rmo/src/rmo.cpp

Lines 213 to 283 in cf85b82

	NumericMatrix sample_cpp(double rate, double rate_killing, double rate_drift, Function rjump, List rjump_arg_list, NumericVector barrier_values) {
	barrier_values = clone(barrier_values);
	if (rate_drift>0.) {
	std::sort(barrier_values.begin(), barrier_values.end());
	} else {
	barrier_values = NumericVector(1, max(barrier_values));
	}
	unsigned int d = barrier_values.size();
	double waiting_time;
	double jump_value;
	double killing_waiting_time;
	double current_value;
	double intermediate_waiting_time;
	double intermediate_value;
	Function do_call("do.call");
	rjump_arg_list.push_back(1, "n");
	std::vector<double> times(1);
	std::vector<double> values(1);
	for (unsigned int i=0; i<d; i++) {
	while (values.back() < barrier_values[i]) {
	waiting_time = ((0. == rate) ? R_PosInf : R::exp_rand()/rate);
	// requires RNGstate synchronisation
	PutRNGstate();
	jump_value = as<double>(do_call(rjump, rjump_arg_list));
	GetRNGstate();
	killing_waiting_time = ((0. == rate_killing) ? R_PosInf : R::exp_rand()/rate_killing);
	if (killing_waiting_time < R_PosInf && killing_waiting_time <= waiting_time) {
	for (unsigned int j=i; j<d; j++) {
	if (rate_drift > 0. && (barrier_values[j] - values.back())/rate_drift <= killing_waiting_time) {
	intermediate_waiting_time = (barrier_values[j] - values.back()) / rate_drift;
	times.push_back(times.back() + intermediate_waiting_time);
	values.push_back(barrier_values[j]);
	killing_waiting_time -= intermediate_waiting_time;
	}
	}
	times.push_back(times.back() + killing_waiting_time);
	values.push_back(R_PosInf);
	} else {
	for (unsigned int j=i; j<d; j++) {
	if (rate_drift > 0. && (barrier_values[j] - values.back())/rate_drift <= waiting_time) {
	intermediate_waiting_time = (barrier_values[j] - values.back())/rate_drift;
	times.push_back(times.back() + intermediate_waiting_time);
	values.push_back(barrier_values[j]);
	waiting_time -= intermediate_waiting_time;
	}
	}
	if (rate > 0.) { // waiting_time < R_PosInf
	times.push_back(times.back() + waiting_time);
	values.push_back(values.back() + waiting_time * rate_drift + jump_value);
	}
	}
	}
	}
	NumericMatrix out(times.size(), 2);
	for (unsigned int i=0; i<times.size(); i++) {
	out(i, 0) = times[i];
	out(i, 1) = values[i];
	}
	colnames(out) = CharacterVector::create("t", "value");
	return out;
	}

Review, simplify, and document the test-implementations in R:

rmo/tests/testthat/helper__lfm_cpp.R

Lines 186 to 203 in cf85b82

	test__rmo_lfm_cpp_R <- function(n, d, rate, rate_killing, rate_drift, rjump_name, rjump_arg_list = list()) { # nolint
	assert_that(is.count(n), is.count(d), is_nonnegative_number(rate),
	is_nonnegative_number(rate_killing), is_nonnegative_number(rate_drift),
	is_positive_number(rate + rate_killing + rate_drift),
	is_rjump_name(rjump_name), is_rjump_arg_list(rjump_name, rjump_arg_list))
	out <- matrix(NA, nrow=n, ncol=d)
	for (k in 1:n) {
	unit_exponentials <- rexp(d)
	cpp_subordinator <- test__sample_cpp_R(rate, rate_killing, rate_drift,
	rjump_name, rjump_arg_list, unit_exponentials)
	out[k, ] <- vapply(1:d, function(x) {
	min(cpp_subordinator[cpp_subordinator[, 2] >= unit_exponentials[[x]], 1])
	}, FUN.VALUE=0.5)
	}
	out
	}

rmo/tests/testthat/helper__lfm_cpp.R

Lines 208 to 258 in cf85b82

	test__sample_cpp_R <- function(rate, rate_killing, rate_drift, rjump, rjump_arg_list, barrier_values) { # nolint
	if (rate_drift>0.) {
	barrier_values <- sort(barrier_values)
	} else {
	barrier_values <- max(barrier_values)
	}
	d <- length(barrier_values)
	times <- 0.
	values <- 0.
	for (i in 1:d) {
	while (last(values) < barrier_values[[i]]) {
	waiting_time <- rexp_if_rate_zero_then_infinity(1, rate)
	jump_value <- do.call(rjump, args=c("n"=1, rjump_arg_list))
	killing_waiting_time <- rexp_if_rate_zero_then_infinity(1, rate_killing)
	if (killing_waiting_time < Inf && killing_waiting_time <= waiting_time) {
	for (j in i:d) {
	if (rate_drift>0. &&
	(barrier_values[[j]] - last(values)) /
	rate_drift<=killing_waiting_time) {
	intermediate_waiting_time <- (barrier_values[[j]] - last(values)) / rate_drift
	times <- c(times, last(times) + intermediate_waiting_time)
	values <- c(values, barrier_values[[j]])
	killing_waiting_time <- killing_waiting_time - intermediate_waiting_time
	}
	}
	times <- c(times, last(times) + killing_waiting_time)
	values <- c(values, Inf)
	} else {
	for (j in i:d) {
	if (rate_drift>0. &&
	(barrier_values[[j]] - last(values))/rate_drift <= waiting_time) {
	intermediate_waiting_time <- (barrier_values[[j]] - last(values)) / rate_drift
	times <- c(times, last(times) + intermediate_waiting_time)
	values <- c(values, barrier_values[[j]])
	waiting_time <- waiting_time - intermediate_waiting_time
	}
	}
	if (rate>0.) { ## waiting_time<Inf # nolint
	times <- c(times, last(times) + waiting_time)
	values <- c(values, last(values) + waiting_time * rate_drift + jump_value)
	}
	}
	}
	}
	cbind("t"=times, "value"=values)
	}

rmo/tests/testthat/helper__lfm_cpp.R

Lines 148 to 160 in cf85b82

	test__rmo_lfm_cpp_independence_R <- function(n, d, rate, rate_killing, rate_drift, rjump_name, rjump_arg_list) { # nolint
	assertthat::assert_that(assertthat::is.count(n), assertthat::is.count(d),
	is_positive_number(rate_drift),
	is_rjump_name(rjump_name), is_rjump_arg_list(rjump_name, rjump_arg_list)) ## dummy test
	out <- matrix(NA, nrow=n, ncol=d)
	for (k in 1:n) {
	unit_exponentials <- rexp(d, 1)
	out[k, ] <- unit_exponentials / rate_drift
	}
	out
	}

rmo/tests/testthat/helper__lfm_cpp.R

Lines 166 to 178 in cf85b82

	test__rmo_lfm_cpp_comonotone_R <- function(n, d, rate, rate_killing, rate_drift, rjump_name, rjump_arg_list) { # nolint
	assertthat::assert_that(assertthat::is.count(n), assertthat::is.count(d),
	is_positive_number(rate_killing),
	is_rjump_name(rjump_name), is_rjump_arg_list(rjump_name, rjump_arg_list)) ## dummy test
	out <- matrix(NA, nrow=n, ncol=d)
	for (k in 1:n) {
	stopifnot(all(rexp(d, 1) >= 0)) ## dummy
	out[k, ] <- rexp(1L, rate_killing)
	}
	out
	}

rmo/tests/testthat/helper__lfm_cpp.R

Lines 6 to 14 in cf85b82

	test__rmo_lfm_cpp_bivariate_R <- function(n, d, rate, rate_killing, rate_drift, rjump_name, rjump_arg_list) { # nolint
	if (rjump_name == "rexp") {
	return(test__rmo_lfm_cpp_bivariate_rexp_R(n, rate, rate_killing,
	rate_drift, rjump_arg_list$rate))
	} else if (rjump_name == "rposval") {
	return(test__rmo_lfm_cpp_bivariate_rposval_R(n, rate, rate_killing,
	rate_drift, rjump_arg_list$value))
	}
	}

rmo/tests/testthat/helper__lfm_cpp.R

Lines 19 to 76 in cf85b82

	test__rmo_lfm_cpp_bivariate_rexp_R <- function(n, rate, rate_killing, rate_drift, jump_rate) { # nolint
	stopifnot(rexp(1L, jump_rate) > 0) ## dummy
	out <- matrix(NA, nrow=n, ncol=2L)
	for (k in 1:n) {
	unit_exponentials <- rexp(2L, 1)
	if (rate_drift>0.) {
	barrier_values <- sort(unit_exponentials)
	} else {
	barrier_values <- max(unit_exponentials)
	}
	d <- length(barrier_values)
	times <- 0.
	values <- 0.
	for (i in 1:d) {
	while (last(values) < barrier_values[[i]]) {
	waiting_time <- rexp_if_rate_zero_then_infinity(1, rate)
	jump_value <- rexp_if_rate_zero_then_infinity(1, jump_rate)
	killing_waiting_time <- rexp_if_rate_zero_then_infinity(1, rate_killing)
	if (killing_waiting_time < Inf && killing_waiting_time <= waiting_time) {
	for (j in i:d) {
	if (rate_drift>0. &&
	(barrier_values[[j]] - last(values))/rate_drift<=killing_waiting_time) {
	intermediate_waiting_time <- (barrier_values[[j]] - last(values)) / rate_drift
	times <- c(times, last(times) + intermediate_waiting_time)
	values <- c(values, barrier_values[[j]])
	killing_waiting_time <- killing_waiting_time - intermediate_waiting_time
	}
	}
	times <- c(times, last(times) + killing_waiting_time)
	values <- c(values, Inf)
	} else {
	for (j in i:d) {
	if (rate_drift>0. && (barrier_values[[j]] - last(values))/rate_drift <= waiting_time) {
	intermediate_waiting_time <- (barrier_values[[j]] - last(values)) / rate_drift
	times <- c(times, last(times) + intermediate_waiting_time)
	values <- c(values, barrier_values[[j]])
	waiting_time <- waiting_time - intermediate_waiting_time
	}
	}
	if (rate>0) { ## waiting_time<Inf # nolint
	times <- c(times, last(times) + waiting_time)
	values <- c(values, last(values) + waiting_time * rate_drift + jump_value)
	}
	}
	}
	}
	cpp_subordinator <- cbind("t"=times, "values"=values)
	out[k, ] <- vapply(1:2L, function(x) min(cpp_subordinator[cpp_subordinator[, 2] >= unit_exponentials[[x]], 1]), FUN.VALUE=0.5) # nolint
	}
	out
	}

rmo/tests/testthat/helper__lfm_cpp.R

Lines 81 to 138 in cf85b82

	test__rmo_lfm_cpp_bivariate_rposval_R <- function(n, rate, rate_killing, rate_drift, jump_value) { # nolint
	out <- matrix(NA, nrow=n, ncol=2L)
	for (k in 1:n) {
	unit_exponentials <- rexp(2L, 1)
	if (rate_drift>0.) {
	barrier_values <- sort(unit_exponentials)
	} else {
	barrier_values <- max(unit_exponentials)
	}
	d <- length(barrier_values)
	times <- 0.
	values <- 0.
	for (i in 1:d) {
	while (last(values) < barrier_values[i]) {
	waiting_time <- rexp_if_rate_zero_then_infinity(1, rate)
	killing_waiting_time <- rexp_if_rate_zero_then_infinity(1, rate_killing)
	if (killing_waiting_time < Inf && killing_waiting_time <= waiting_time) {
	for (j in i:d) {
	if (rate_drift>0. &&
	(barrier_values[[j]] - last(values))/rate_drift<=killing_waiting_time) {
	intermediate_waiting_time <- (barrier_values[[j]] - last(values)) / rate_drift
	intermediate_value <- barrier_values[[j]]
	times <- c(times, last(times) + intermediate_waiting_time)
	values <- c(values, intermediate_value)
	killing_waiting_time <- killing_waiting_time - intermediate_waiting_time
	}
	}
	times <- c(times, last(times) + killing_waiting_time)
	values <- c(values, Inf)
	} else {
	for (j in i:d) {
	if (rate_drift>0. &&
	(barrier_values[[j]] - last(values))/rate_drift <= waiting_time) {
	intermediate_waiting_time <- (barrier_values[[j]] - last(values)) / rate_drift
	intermediate_value <- barrier_values[[j]]
	times <- c(times, last(times) + intermediate_waiting_time)
	values <- c(values, intermediate_value)
	waiting_time <- waiting_time - intermediate_waiting_time
	}
	}
	if (rate>0.) { ## waiting_time<Inf # nolint
	times <- c(times, last(times) + waiting_time)
	values <- c(values, last(values) + waiting_time * rate_drift + jump_value)
	}
	}
	}
	}
	cpp_subordinator <- cbind("t"=times, "values"=values)
	out[k, ] <- vapply(1:2L, function(x) min(cpp_subordinator[cpp_subordinator[, 2] >= unit_exponentials[[x]], 1]), FUN.VALUE=0.5) # nolint
	}
	out
	}