When I simulate normal data in R, I make sure that the sample have the exact mean and sd of the sampling distribution:
x = scale(rnorm(n))*sd + mean.
I want to do the same for binomial data, making the sample express the near-exact probability that they were generated from. Of course it can't be exact when the probability is continuous and the sample is discrete but something that gets pretty close would be nice
x = rbinom(n, 18, 0.5) can potentially give samples where an MLE estimate would indicate a probability of p=0.2 or p=0.8 which is pretty far from p=0.5.
Purpose: I'm building a Bayesian model where I infer a binomial rate from a small sample. To test that the model works, I'd like to simulate well-specified data, in order to diagnose whether a strange inferential result is due to chance in the simulation or in the model.