Skip to content

Instantly share code, notes, and snippets.

@dickoa
Last active August 27, 2023 16:31
Show Gist options
  • Save dickoa/eed351acfcdf7383ee900491efcd1ff6 to your computer and use it in GitHub Desktop.
Save dickoa/eed351acfcdf7383ee900491efcd1ff6 to your computer and use it in GitHub Desktop.
Brewer method unequal probabilities speed up
#include <R.h>
#include <Rmath.h>
void C_up_brewer(double* pik, int* N, double* sb) {
double n = 0.0;
double a;
double u;
int j;
// Initialize sb and calculate sum of pik (n)
for (int i = 0; i < *N; i++) {
sb[i] = 0;
n += pik[i];
}
// Random seed
GetRNGstate();
// Main loop
for (int i = 1; i <= (int)n; i++) {
double p_sum = 0.0;
double* p = (double*) malloc(*N * sizeof(double));
// Calculate a
a = 0.0;
for (int k = 0; k < *N; k++) {
a += pik[k] * sb[k];
}
// Calculate p
for (int k = 0; k < *N; k++) {
p[k] = (1 - sb[k]) * pik[k] * ((n - a) - pik[k]) / ((n - a) - pik[k] * (n - i + 1));
p_sum += p[k];
}
// Generate a random number from a uniform distribution
u = runif(0, 1);
// Find the first index j such that p[j] - u > 0
double cumsum = 0.0;
for (j = 0; j < *N; j++) {
cumsum += p[j]/p_sum;
if (cumsum - u > 0) {
break;
}
}
// Update sb
sb[j] = 1;
// Free the dynamically allocated memory for p
free(p);
}
PutRNGstate();
}
up_brewer <- function(pik, eps = 1e-06) {
if (any(is.na(pik)))
stop("there are missing values in the pik vector")
list <- pik > eps & pik < 1 - eps
pikb <- pik[list]
s <- pik
N <- as.integer(length(pikb))
sb <- vector(mode = "double", length = N)
r <- .C("C_up_brewer", pikb, N, result = sb)
s[list] <- r$result
s
}
@dickoa
Copy link
Author

dickoa commented Aug 27, 2023

In linux, you can compile up_brewer.c from R using

system("rm up_brewer.o up_brewer.so")
system("R CMD SHLIB -lRmath up_brewer.c")

And don't forget to load .so file before running up_brewer function in R.
dyn.load("up_brewer.so")

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment