Calculate the probability of observing n elements in m draws from an uneven frequency distribution, in order but consecutive or not I want to calculate the probability to observe a certain sequence of n elements in order, but not necessarily consecutively, after drawing m times from an uneven frequency distribution.
For example:
Frequency distribution:
A- 0.5 B- 0.4 C- 0.05 D-0.05
What's the probability to observe the sequence 'AAA' if you draw 10 times from this distribution? A good draw could have consecutive A's or not (AAACCCCCCC or ACCCCACCCA are both valid). What about ACB?
The goal is to filter out matches that are too common and short from sequence alignment. For example, aligning AAA with ABCABCACBCBC does not mean much as the A's are frequent in the distribution.
One extra problem is that simulation approaches are not possible here, as I have multiple millions of such comparisons to score.
 A: You can write a dynamic program to compute these values in time proportional to the length of the probe string and the number of draws.
For a probe string $s=s_1s_2\cdots s_k$ where letter $s_i$ has probability $p_i$, let $f(s,m)$ be the desired probability. Drawing from the distribution, considering the chance that the first letter in the target matches $s_1$, and writing $s^\prime=s_2s_3\cdots s_k$ for the suffix, it is immediate that $$f(s,m)=p_1f(s^\prime,m-1)+(1-p_1)f(s,m-1).$$ For initial conditions use $f(s,0)=0$ for nonempty $s$ and $f(s, n)=1$ for empty $s$.
The calculations are rapid.
Examples: $f(AAA,10) = 121/128$ and $f(ABC,10)=1247150641/5120000000.$

Here is working R code to illustrate.
f <- function(p, n) {
  p <- rev(p)
  k <- length(p)
  a <- rbind(1, matrix(0, k, n+1))
  for (i in 1:k) {
    for (j in i:n+1) 
      a[i+1, j] <- p[i] * a[i, j-1] + (1-p[i]) * a[i+1, j-1]
  }
  return(a[k+1, n+1])
}

For example,
a <- .5; b <- 0.4; c <- 0.05; d <- 0.05
f(c(a,a,a), 10)


[1] 0.9453125

print(f(c(a,c,b), 10), digits=16)


[1]  0.2435841095703125


For symbolic solutions, here is working Mathematica code:
f[x_, n_] := f[x, n] = Expand[x [[1]] f[Rest@x, n - 1] + (1 - x [[1]]) f[x, n - 1]];
f[x_, 0] /; Length[x] > 0 := 0;
f[{}, n_] /; n >= 0 = 1;

For example,
f[{a,b,c}, 5]


$a^3 b c+a^2 b^2 c+a^2 b c^2-5 a^2 b c+a b^3 c+a b^2 c^2-5 a b^2 c+a b c^3-5 a b c^2+10 a b c$

f[{a,b,c}, 10] //Length


120

(The general answer is a polynomial with 120 terms.)
f[{10, 8, 1}/20, 10]


$ \frac{1247150641}{5120000000}$

N[f[{10, 8, 1}/20, 10], 16]


$0.2435841095703125$

(This agrees with the R solution above.)
