Case scenario: a retro-virus infects a healthy cell. The virus programs the cell to brew little viruses, at a rate of 0.5 per-sec, until finally the cell bursts when the number of virus inside it is 5. How to model this?

In Binomial, the random variable represents the number of successful trials obtained when throwing a coin a certain number of trials, at a certain probability of success per trial.

I want a distribution whose random variable is the number of trials (coin tosses) that were necessary to perform, given a certain number of successful trials and a certain probability per trial.

I am not even sure how I would write down the probability mass function.

Is there such a distribution? Nothing rings a bell here: https://en.wikipedia.org/wiki/Category:Discrete_distributions

There are related questions to this one --- such as this one: How many trials until I each my desired outcome --- but no-one mentioned a distribution, or if any exists.

Just to make it clear, the random generator for a random variable of such a distribution would look like this in R:

rmy <- function(s, p) {
    i <- n <- 0
    while(i != s) {
       i <- i+rbinom(1,1,p)
       n <- n+1

Thank you ! ps: sorry if the text was a little flowery, but it helps me think, since I am a junior mathematician hehe.

  • $\begingroup$ Let me see if I understand the process. Each second I flip a (fair) coin to see if I produce a virus. Since I have a $1/2$ probability of producing a single new virus, the average virus production rate is $1/2$. I stop when I produce five viruses. And you want to know the expected number of coin flips, is that right? $\endgroup$ – Brian Tung Apr 14 '15 at 22:23
  • $\begingroup$ Have you looked at Poisson Distribution? $\endgroup$ – rightskewed Apr 14 '15 at 22:46
  • $\begingroup$ @BrianTung, I want to know if there is a theoretical distribution. Expected value is easy, it is E=m/p=5/0.5=10, in this case, where m is number of successful trials and p is probability of success. $\endgroup$ – rpmcruz Apr 15 '15 at 8:05
  • $\begingroup$ @rightskewed, hmm what do you have in mind exactly? Using $\lambda=m/p$; testing in R rpois(10, 5/0.5) shows it sampling values below <5, which is impossible in this model. $\endgroup$ – rpmcruz Apr 15 '15 at 8:08
  • $\begingroup$ Just to make it clear: I want to know if there is already a popular theoretical distribution for what I want. If anyone works to work out the probability mass function that would be superb too. :P Might make sense to migrate it to: stats.stackexchange.com . I did not remember about the statistics forum ... $\endgroup$ – rpmcruz Apr 15 '15 at 8:12

I figure this one out. :)

I can model it using a Negative Binomial: https://en.wikipedia.org/wiki/Negative_binomial_distribution

First, let us change the values of my case scenario, just to make it clearer. "Case scenario: a retro-virus infects a healthy cell. The virus programs the cell to brew little viruses, at a rate of 0.2 per-sec, until finally the cell bursts when the number of virus inside it is 5. How to model this?"

We can model number of failures $Y$ as $Y\sim\mathcal{NB}(5,0.2)$. That answers the question, how many failed trails do we have, when we need 5 successful at a probability rate of 0.2. But we do not want failed trials, we want total trials, and total trials = failed trials + successful trials. We know successful trials, which is 5, so our random variable $X$ is such that $X\sim5+\mathcal{NB}(5,0.2)$.

In fact, comparing the random generator function I proposed in the question with the negative binomial random generator (with this adjustment):

hist(sapply(1:1e5, function(x) rmy(5, 0.2)))
hist(5+rnbinom(1e5, 5, 0.2))

testing distribution

All functions mean, sd and summary are consistent as well.


My 2 cents here :

Using a Negative Binomial distribution can be a very good approximation. However :

1) You totally remove the case r=0 (No failures).

2) Also, you're not answering your initial question. You're computing the distribution of the #trials before achieving the known #fails (or #success). It's different than #trials knowing the #fails (or #success), as you can have fails following the last success !

I believe the right answer involves Bayesian probabilities :

  • We consider that N is a random variable describing the # of trials (# of Bernoulli trials). Let's assume that N is following an Exponential Distribution $ N \thicksim Exp(\lambda) $
  • We then consider that X is a random variable describing the number of successes after N=n trials, following a Binomial distribution $ N \thicksim Binom(p, N) $

Let's compute $ P(n|k) $ thanks to the Bayes formula :

$$ \frac{1}{P(N=n|X=k)} = \frac{P(X=k)}{P(X=k|N=n)P(N=n)} $$

And by writing $ P(k) = \sum_{m=k}^{\infty}{P(X=k|N=m)P(N=m)} $

$$ = \sum_{m=k}^{\infty}\frac{P(X=k|N=m)P(N=m)}{P(X=k|N=n)P(N=n)} $$

And because $P(X=k|N=n)$ and $P(X=k|N=m)$ are just Binomial distributions :

$$ = \sum_{m=k}^{\infty}\frac{m!k!(n-k)!}{n!k!(m-k)!} \frac{p^k}{p^k} \frac{(1-p)^{m-k}}{(1-p)^{n-k}} \frac{e^{-\lambda.m}}{e^{-\lambda.n}}$$

$$ = \sum_{m=k}^{\infty}{[e^{-\lambda}(1-p)]^{m-n}}$$

And by re-indexing the Serie with $l=m-k <=> m=l+k$

$$ = \sum_{l=0}^{\infty}{[e^{-\lambda}(1-p)]^{l+k-n}}$$

$$ \frac{1}{P(N=n|X=k)} = \frac{[e^{-\lambda}(1-p)]^{k-n}}{1 - e^{-\lambda}(1-p)}$$

$$ P(N=n|X=k) = \frac{1 - e^{-\lambda}(1-p)}{[e^{-\lambda}(1-p)]^{k-n}}$$

And now you have a probability $P>0$ for $k=0$ ;-)

Edit: The choice of the prior (Exponential) is totally arbitrary and leads to a lot of simplifications. If anyone wants to attempt the calculation with another prior, feel free :-)


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.