Is there a more efficient method for creating a probability distribution?

Question

I have written a short function in R to estimate the expected number of mutations that will be observed in a set of DNA sequences. The parameters are the mutation rate (x), the length of the DNA sequence (y) and the number of stretches of DNA being observed (z).

I want to simulate how many mutations one expects to observe in real data, given these underlying parameter values, so I do random sampling using these values and repeat this 100 times to get a distribution of the observed number of mutations.

In the example below the mutation rate is 1x10-8, in a 100000 nucleotide DNA sequence, observed in 1000 individuals:

x = 1/100000000
y = 100000
z = 1000

mut_counter <- function(x,y,z) { 
  expect_muts <- replicate(100, {
    sum(observed_muts <- sample(0:1, (y*z), prob=c(1-x,x), replace = TRUE)>0)
  })
  return(expect_muts)
}

The output is a list with the number of observed mutations in each of the 100 replicates. Here are the first 3 rows (only 1 mutation is observed in one replicate):

 [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [21] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
 [41] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

This random sampling is slow when the parameter values become very large. Is there a way to achieve the same outcome using a more efficient method? For example I was wondering if sampling from a poisson distribution could be used instead or if that is not suitable for this case.

Thank you for you assistance. Please let me know if I can be more clear.

Answer 1

You could use the exact binomial (and it will likely be considerably faster), but the Poisson will be highly accurate.

Indeed, this is the ideal situation for the Poisson.

The expected number is simply x*y*z (i.e. 1 for your numbers). Your expected_muts values should be essentially Poisson(x*y*z) (i.e. Poisson(1)); if you need a sample from those it would be vastly faster to simulate them directly. On the other hand if you just need to know the probabilities from the distribution, rather than sample 100 times you can produce the probabilities of each outcome out as far as you need:

[In the individual sample (observed_muts) the expected number is $10^{-8}$; if you really needed those for anything you could do those as either Poisson($10^{-8}$) or as binomial; it should make no appreciable difference.]

It's not clear exactly what you're trying to do with these numbers, which makes it harder to give more suitable advice.

---

R code that gives a similar image to the above

x = 1/100000000
y = 100000
z = 1000
xi = 0:6
pb=dbinom(xi,10^8,10^-8)
pp=dpois(xi,1)
plot(pb~I(xi-.05),type="h",xlab="x",ylab="p")
points(pp~I(xi+.05),col=2,type="h")