Using R's internal sample function in CSnippet

[trh3]

Dr. King,

I'm hoping that this is a fairly simple question to answer, but I'm trying to set up some ordinal measurement state space models in pomp, and am running into trouble figuring out how to integrate R's internal sample function. More specifically, I really don't know which of the functions related to sample to use. I'm fairly new to using R's internal C functions, so I'm likely missing something quite obvious. Happy to provide the code that I am working on if that helps clarify the issue.

Thank you!
Teague

[kingaa]

Have you looked at the documentation on the C API? In particular, do you know that there is a C entry point for sample? If you don't find the answer there, do include the code and I can try to help you figure it out.

[trh3]

Aha, I was looking at the R source code and trying to pull functions from there, rather than looking at the actual documentation, thanks for pointing me to that! That said, there doesn't seem to be an entry point for sample, none described in the documentation at least (though again, I could be missing something). I did manage to use rmultinom and then check the sampled vector, which is working well for my purposes. Relevant C code for others to reference below:

/*Hardcoded "sample" for 6 categories*/
/* Allocate final integer object 
int y1_temp = 0;
/* Allocate the vector that will contain the sample */
int perm1[6] = {0,0,0,0,0,0};
/* Sample a single draw from multinomial distribution (p1diff calculated previously) */
rmultinom(1, &p1diff[0], 6, &perm1[0]);
/* Transform from vector of 0/1 to integer (base 1) index */
for(int i = 0; i < 6; i++){
  if(perm1[i]== 1){
    y1_temp = i+1;
  }
}

I'm sure there is a bit of a simpler/tidier way, but this gets me what I am looking for!

[kingaa]

Yes, this will work for sampling with replacement, but I believe it may not be very efficient. Looking into the innards of the rmultinom function, you see that it requires repeated calls to a binomial sampler, each of which requires a minimum of one runif. In principle, to draw a single random sample from a set of choices should require no more than 1 runif call. See, for example, the source code beginning at line 312 of random.c.

[trh3]

I thought the multinomial sampler might not be very efficient! The structure of this model makes it amenable to the runif call/ check against cumulative distributions. I'll modify and post that code here for reference.

[trh3]

If p1 is an array of cumulative probabilities ordered descending with a 0 appended (in this case, 6 categories, array length of 7), then one can sample from the set of choices as follows:

double draw1 =  unif_rand();
int y1_temp = 0;
for(int i = 6; i > 1; i--){
  if(draw1 > p1[i]){
    y1_temp = i;
  }
}

[kingaa]

See also the code for nosort_resamp in pomp.

[trh3]

Fantastic thank you! This has been very helpful!