3
$\begingroup$

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.

Improve this question
$\endgroup$
10
  • 3
    $\begingroup$ Hmm. That's kinda a strange request but .. Binomial gives integers so with fixed n, you can't always find a set of integers that will exactly hit a mean, say Pi. Now, if that mean is from another binomial than you can obviously but perhaps not in a different way. My guess is that you may not have to do this. What's the purpose? $\endgroup$
    – Dave31415
    Mar 24, 2014 at 15:22
  • 5
    $\begingroup$ This question still seems under-specified. Why, for instance, would it not suffice to generate c(rep(1, round(n*p)), rep(0, n-round(n*p)))? That would seem to satisfy all the requirements you have set out. $\endgroup$
    – whuber ♦
    Mar 24, 2014 at 15:53
  • 7
    $\begingroup$ "When I simulate normal data in R, I make sure that the sample have the exact mean and sd of the sampling distribution" --- Except that by doing so, you guarantee it's not actually normally distributed. Is that actually what you want here? And if you do want to absolutely guarantee it's not actually a random sample from a binomial, why bring up 'binomial' at all? $\endgroup$
    – Glen_b
    Mar 24, 2014 at 22:40
  • 1
    $\begingroup$ (ctd) Though I guess the easy way to do that would be to multiply the pmf by $n$ and round, and then add or remove a few values that leave exactly $n$, by adding to values that are most under-represented or subtracting from ones that are over, in such a way as to get the desired moments close to right. $\endgroup$
    – Glen_b
    Mar 24, 2014 at 23:25
  • 4
    $\begingroup$ Jonas, the deterministic sample I specified is as close to Binomial as you possibly can get for any given $n$, so it is not evident what you would mean by "not the best approximation." (In fact, my example is a very simple illustration of the general procedure proposed by @Glen_b.) It may be worth bluntly stating that what you are asking for looks like a poor way of carrying out your intended investigation. Instead, generate random samples according to simpler procedures, without any constraints, and analyze the effects of sampling variability on your simulation. $\endgroup$
    – whuber ♦
    May 26, 2014 at 15:00

2 Answers 2

Reset to default
2
$\begingroup$

I'm not sure whether this is actually advisable, but it should be straightforward to generate. What you are asking for, essentially, is an underdispersed binomial distribution. You can get this by sampling (with replacement, if you want more than 1 value) from a vector of the integers 0:size, where you specify a set of underdispersed probabilities. You just have to figure out what the probabilities are that you want. First, consider this figure:

enter image description here

From this you can see that the probabilities of each possible value from a binomial distribution can be matched by a normal distribution with mean $n\pi$ and variance $n\pi(1-\pi)$. Thus, you can make an underdispersed version by using appropriately scaled densities from a normal distribution with the same mean but a smaller SD. Imagine that you want the SD to be cut in half, then: 

set.seed(1773)
hSD.norm  = dnorm(0:18, mean=9, sd=sqrt(18*.25)*.5)
ud.probs  = hSD.norm/sum(hSD.norm)
N         = 10000
vals      = sample(0:18, size=N, replace=TRUE, prob=ud.probs)

enter image description here

Improve this answer
$\endgroup$
0
1
$\begingroup$

One particularly blunt way of doing it is to keep generating samples until you get a sample that you like (i.e. sample mean close to population mean).

N <- 100
n <- 18
p <- 0.5
e <- 0.01

max.attempts <- 10000

candidate <- rbinom(N, n, p)
attempt <- 1

while (abs(mean(candidate) - n*p) >= e &&
        attempt <= max.attempts) {
  candidate <- rbinom(N, n, p)
  attempt <- attempt + 1
}

hist(candidate,
     main=paste("Binomial sample with mean", mean(candidate)))

It is straightforward to extend this to meet a requirement on the standard deviation or any other moments. Just make sure that it's actually possible to get as close as you specify to the given mean.

Improve this answer
$\endgroup$
1
  • 3
    $\begingroup$ This is a reasonable (and sometimes effective) approach. However, it carries two risks that would be worth discussing. First, it can take a huge number of trials unless the threshold e is made adequately large. Second, the resulting sample might look nothing at all like a binomial distribution! A modification of your approach that uses the methods suggested in comments by @Glen_b would resolve the second issue. $\endgroup$
    – whuber ♦
    Apr 24, 2014 at 18:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.