1
$\begingroup$

Suppose I am generating 10000 data points with probability P out of N total trials, but I know that P x N can never be below M, where M < N. How would I go about restricting my simulation to never have P x N < M?

Specifically, P=.9 and N=30 and M = 22. If I simulate 10000 values, I will almost inevitably have some values less than 22. My current method is to take the values that were less than 22 and re-populate them using the same binomial values, and repeating that until I have no values less than 22. My current R code is here

x=rbinom(10000, 30, .9) # yields some values less than 22.
l=length(x[x<22]) # Values less than 22 (a handful)
while({l}) { x[x<22] <- rbinom(l,30,.9); l=length(x[x<22]) }

Also, is this applicable to other distributions? Particularly the Beta distribution when probabilities are in mind.

$\endgroup$
5
  • 1
    $\begingroup$ Do you care that the expected value of your truncated distribution is greater than $NP$? (in your example it would be $27.0127$ rather than $27$, but the gap would change if you used a cut-off higher than $22$) $\endgroup$
    – Henry
    Commented Apr 6, 2020 at 17:29
  • $\begingroup$ Code like x <- rbinom(20000, 30, 0.9); x <- x[x >= 22][1:10000] might be faster $\endgroup$
    – Henry
    Commented Apr 6, 2020 at 17:33
  • 3
    $\begingroup$ Your problem is difficult to decipher, because neither $P$ nor $N$ depend on the data: you specify them as part of your program. The evident answer is not to specify values for which $PN\lt M.$ What more is there to say?? Do you mean that you do not want any of the values to be less than $M$? That would be a truncated Binomial distribution. If so, your method is a special case of "rejection sampling," which you can investigate here on CV with a search. $\endgroup$
    – whuber
    Commented Apr 6, 2020 at 18:38
  • $\begingroup$ @whuber My example was notional, but this is testing data. M is actually the current successful testing rate, and we are adding N-M additional tests to the data. So the reason for the restriction is to limit the distribution to be at least the minimum number of successes. $\endgroup$
    – Sumner18
    Commented Apr 6, 2020 at 22:20
  • $\begingroup$ @whuber See my other question at stats.stackexchange.com/questions/458331/… where this code is coming from. $\endgroup$
    – Sumner18
    Commented Apr 6, 2020 at 22:32

1 Answer 1

1
$\begingroup$

You seek $X \sim \mathsf{Binom}(n = 30, p = .9)$ truncated to avoid values smaller than 22.

Without the truncation, you could use the inverse CDF (quantile) method to get $Y \sim \mathsf{Binom}(n = 30, p = .9)$ in R as follows;

y = qbinom(runif(10^5, 0, 1), 30, .9)

In order to avoid the unwanted outcomes, observe that $P(Y \le 21) = 0.007783619.$ So truncate the standard uniform distribution to $\mathsf{Unif}(0.007783619, 1)$ in order to avoid generating the unwanted observations.

pbinom(21, 30, .9)
[1] 0.002019829

Suggest R code:

set.seed(2020)
x = qbinom(runif(10^5, pbinom(21,30,.9) ,1), 30, .9)
table(x)/10^5
x
     22      23      24      25      26      27      28      29      30 
0.00587 0.01808 0.04787 0.10307 0.17904 0.23915 0.22421 0.14114 0.04157     
pdft= dbinom(22:30, 30, .9)
pdf = pdft/sum(pdft);  round(pdf, 5)
[1] 0.00578 0.01808 0.04746 0.10251 0.17742 0.23657 0.22812 0.14159 0.04248

hist(x, prob=T, br=seq(21.5, 30.5, by=1), col="skyblue2")
points(22:30, pdf, pch=10, col="red") 

enter image description here

Note: The method seems to work for other truncated distributions as well. It is simple and efficient. I have not seen it before, but I would not be surprised to hear that it is in general use.

$\endgroup$
5
  • $\begingroup$ Does this address the concern from @Henry in the comments above that the Expected Value would change and be greater than NP with different restrictions? $\endgroup$
    – Sumner18
    Commented Apr 6, 2020 at 22:27
  • $\begingroup$ Inevitably a one-sided truncation will change the mean and median. A half-normal distribution certainly does not have mean 0. I'd say @Henry is pointing out an obvious fact, not voicing a 'concern'. [Indeed, R code 'sum((22:30)*pdf)' does return 27.01271.] $\endgroup$
    – BruceET
    Commented Apr 6, 2020 at 22:33
  • $\begingroup$ Is there some way to weight the values to maintain a constant mean and median across changing restrictions? $\endgroup$
    – Sumner18
    Commented Apr 6, 2020 at 22:36
  • $\begingroup$ Seems self-deceptive. Why would you want to do that? $\endgroup$
    – BruceET
    Commented Apr 6, 2020 at 22:39
  • 1
    $\begingroup$ I don't think that it would be important. I'm essentially using this truncated distribution to get possible success rates of additional tests to an already guaranteed number of successes (the restriction). I'm struggling to see whether or not these weights would be important to the distribution. $\endgroup$
    – Sumner18
    Commented Apr 6, 2020 at 22:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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