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

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):

  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
 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.

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")

• Thanks very much for your answer. The aim is to get an expectation of how many mutations we can expect to observe in a stretch of DNA if we sequence x number of individuals. This can help to determine how much sequencing to do if we want to observe at least x mutations to have power for subsequent tests of genome evolution. – user964689 Mar 22 '17 at 13:21
• Would you share the code to make the plot above? Now I am sampling from the Poisson distribution to get a distribution of the values we might observe from real data but it would be nice to also present the probability distribution as you have. My current code: rpois(100,(xyz)), which I then plot as a histogram. – user964689 Mar 22 '17 at 19:47
• See edit... ... 100 observations isn't nearly enough for what I'd call an accurate idea, but you should just draw the probability function – Glen_b Mar 23 '17 at 1:00
• Thanks very much. Yes I was just starting with 100 to see that the code works before scaling up. Seems to work well with the Poisson approach. – user964689 Mar 23 '17 at 16:48
• Note that my graph (using dbinom and dpois) is what you should see as $n\to\infty$ with rbinom or rpois – Glen_b Mar 23 '17 at 19:48