I have a feeling there is a straight forward answer to my question, but I'm not sure what the appropriate term to search for is. I'm trying to get a sense of what the probability is of finding some outcome as a subset of a random sample. The challenge as I see it is that the possible subsets of the total sample are not independent.

As a simple example, imagine I flipped a coin 100 times and recorded the results. I'm trying to ask a question like "what is the probability that somewhere in these results there is a run of 5 heads in row (or some other specific sequence)". The probability of flipping a coin 5 times and getting heads each time is pretty easy to calculate, and you have numerous opportunities within 100 total flips, so I can imagine using the binomial distribution, but the question would be exactly how many 'independent' opportunities you have within all the flips.

Naively, you might say that there are 95 chances - you can get all heads over the first five flips, or failing that, over the flips numbered 2 through 6, or 3 through 7, etc - but that's clearly not right, because these subsets of flips overlap each other. If you know that the 5th flip was tails, then you know that all 5 runs of flips that contain it have been 'spoiled'.

On the other end, you might imagine segmenting the 100 flips into 20 completely independent (i.e non-overlapping) runs of flips, but this would seem to undercount your probability - the overlapping runs clearly add some opportunity of finding the sequence of heads.

The actual problem I'm dealing with that caused me to have this question is only very slightly more complicated than the dummy example with the coin flips, because I'm looking in two dimensions (looking for a particular grid of values in an array of random variables). Somehow it feels like there ought to be a way to account for the overlapping samples and come to an intuition for how often you expect to find some specific sequence in the larger sample.

  • 1
    $\begingroup$ The one-dimensional problem has already been tackled on this website in several ways, for instance Probability of a similar sub-sequence of length X in two sequences of length Y and Z and A fair die is rolled 1,000 times. What is the probability of rolling the same number 5 times in a row?. In order to tackle your 2d problem one can use similar approaches, but it would help to explain it in more detail. $\endgroup$ May 2 at 12:40
  • $\begingroup$ You do not have independence. And you have $96$ rather than $95$ possible strings of $5$, so the expected number of appearances including overlaps of any given sequence of five is $96 \times 2^{-5}=3$, while the expected number excluding overlaps depends on the particular given sequence: for $5$ heads in a sequence on $n$ flips it seems to be about $\frac{1}{31}(\frac n2 - \frac{49}{31})$ so about $1.56$ when $n=100$. This does not help with the probability of ever seeing $5$ heads, which is empirically about $1-0.945595\times 0.982974^n$ for large $n$ about $0.81$ when $n=100$ $\endgroup$
    – Henry
    May 2 at 12:57
  • $\begingroup$ @SextusEmpiricus Thanks for those links, very helpful. Reading up on the Markov Chains now. My 2d problem is basically the same sort of thing, except it's a table of what is essentially noise, and the question is, what is probability that some sub-grid of it has some pattern. So, in my, say, 100 by 100 table of random 1s and 0s, what's the probability I can find, say, a 5 by 5 grid of just 1s. $\endgroup$
    – Max Miller
    May 2 at 13:39
  • $\begingroup$ @MaxMiller Markov Chains are actually the thing where I have doubts whether these can help you (since they are more one dimensional). $\endgroup$ May 2 at 13:41

1 Answer 1


I'm not completely satisfied with this answer, but I think it's better than no answer. I'll attempt to answer the question posed in the comments: in a 100-by-100 matrix of independent Bernoulli trials ($M$), what is the probability of finding a 5-by-5 submatrix, $X$, composed entirely of successes? I use a 0.5 probability of the Bernoulli trials returning 1, but the methodology should work for any probability. The code below is in R.

First, I used a couple Rcpp functions. The first is a fast Rcpp function saved in "C:/temp/findSquare.cpp". It counts how many matrices contain a square submatrix of size n composed entirely of 1s. Each column of the input matrix x holds the values of an nr-by-nc matrix.

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
int findSquare(IntegerMatrix& x, const int& nr, const int& nc, const int& n) {
  const int iter = x.nrow();
  int out = 0;
  for (int k = 0; k < iter; k++) {
    int jj = nc;
    int cnt = 0;
    int found = 0;
    for (int j = 0; j < nc; j++) {
      cnt = (cnt + 1)*x(k, j);
      x(k, j) = (cnt >= n) + 1;
    for (int i = 1; i < nr; i++) {
      int cnt = 0;
      for (int j = 0; j < nc; j++) {
        cnt = cnt*x(k, jj) + x(k, jj);
        x(k, j) = x(k, j)*(cnt >= n) + 1;
        found += x(k, j) > n;
    out += found > 0;

The second is for matrix multiplication since matrix multiplication in base R is so slow on my setup. It is saved in "C:/temp/eigenMult.cpp"

// [[Rcpp::depends(RcppArmadillo, RcppEigen)]]
#include <RcppArmadillo.h>
#include <RcppEigen.h>

// [[Rcpp::export]]
SEXP eigenMapMatSquare(const Eigen::Map<Eigen::MatrixXd> X){
  Eigen::MatrixXd X2 = X * X;
  return Rcpp::wrap(X2);

// [[Rcpp::export]]
SEXP eigenMapMatMult(const Eigen::Map<Eigen::MatrixXd> A, Eigen::Map<Eigen::MatrixXd> B){
  Eigen::MatrixXd C = A * B;
  return Rcpp::wrap(C);

Finally, an R function for raising a matrix to a power using "eigenMult.cpp":

matpow <- function(x, pow) {
  if (pow == 0L) {
  } else if(pow %% 2L) {
    nextpow <- (pow - 1L)/2L
    if (nextpow) {
      eigenMapMatMult(x, Recall(eigenMapMatSquare(x), nextpow))
    } else x
  } else {
    Recall(eigenMapMatSquare(x), pow/2L)

We can get a credible interval via Monte Carlo simulation:

# Monte Carlo - take 100 sets of 10k samples
system.time(counts <- sum(replicate(100, findSquare(matrix(sample(0:1, 1e8, TRUE), 1e4), 100L, 100L, 5L))))
#>   user  system elapsed
#> 641.66   27.36  669.52
# 90% credible interval
(pSimBounds <- qbeta(c(0.05, 0.95), sum(counts) + 0.5, 1e6 - sum(counts) + 0.5))
#> [1] 0.0002212588 0.0002728769

These bounds are fairly tight, but they took more than 11 minutes to compute. A simple Markov chain will give quick and easy bounds. Instead of a 100-by-100 matrix, look for $X$ in a 100-by-5 matrix. The probability of getting 1 for all five values in a row is $2^{-5}$. The transition matrix is then:

m <- cbind(matrix(c(rep(1 - 1/32, 5), 0), 6), diag(1/32, 6)[,-6])
m[6, 6] <- 1

Get the probability of finding $X$:

(p1 <- matpow(m, 100)[1, 6])
#> [1] 2.772544e-06

The full matrix has 20 independent sets of 5 columns. A lower bound could be the probability of finding $X$ only in those 20 independent sets. For an upper bound, there are 96 sets of 5 columns in the full matrix, but they aren't independent: if the $X$ is not in columns 1 through 5, it is less likely to be in columns 2 through 6, etc. Treating the 96 sets as independent will over-estimate the probability of finding $X$ in $M$:

(pLower <- 1 - (1 - p1)^20)
#> [1] 5.544942e-05
(pUpper <- 1 - (1 - p1)^96)
#> [1] 0.0002661292

The lower bound is not very impressive, but the upper bound is tighter than what was calculated from 1M samples.

For the final estimate, let's examine what happens if we start with 5 columns and repeatedly add one more column. For six columns, the Markov chain needs to track what is happening in columns 1:5 as well as 2:6. Each set can be in state 0 through 4 (the number of consecutive rows of 5 1s), resulting in 25 states plus 1 absorbing state. As an example, the sequence 0,1,1,1,1,1 in a row will transition columns 1:5 to state 0 and increment the state of columns 2:6 by 1 (unless the system is in the absorbing state).

If we start with 5 columns and add 1 column at a time, we can calculate the probability of finding $X$ in the last 5 columns conditioned on $X$ not existing in the matrix without the last column. As more columns are added, this probability quickly converges, which allows us to estimate the probability of finding $X$ in $M$. The size of the transition matrix grows quickly as more columns are added, so I stop at 11.

# more involved estimate from multiple Markov chains
library(Rfast) # for eachrow and rowMaxs functions

fTrans <- function(dt, mTrans, transk, n, nn) {
  # iteratively build out state changes
  nmTo <- paste0("to", 1:nn)
  nmFrom <- colnames(dt)
  nmAll <- c(nmFrom, nmTo)
  dtStates <- setnames(as.data.table(rep(list(integer(0)), nn*2L + 1L)), c(nmAll, "k"))
  while (nrow(dt)) {
    ldt <- vector("list", nrow(dt) + 1L)
    ldt[1:nrow(dt)] <- lapply(
      function(i) {
        m2 <- eachrow(mTrans, as.numeric(dt[i])) + mTrans
        m2[rowMaxs(m2, TRUE) == n, ] <- n
          setnames(dt[rep(i, nrow(mTrans))], nmFrom),
          setnames(as.data.table(m2), nmTo)
        )[, k := transk][, .(k = sum(k)), nmAll]
    ldt[[nrow(dt) + 1L]] <- dtStates
    dtStates <- rbindlist(ldt)
    dt <- rbindlist(
        unique(dtStates[, 1:nn])[, k := 2L],
        unique(dtStates[, (nn + 1L):(2L*nn)])[, k := 1L]
      ), FALSE
      , .(k = sum(k)), nmFrom
      k == 1L & from1 < n
      , k := NULL

pGridInArray <- function(n, dims, p = 0.5) {
  # main function
  if (dims[1] < dims[2]) dims <- dims[2:1]
  n1 <- n - 1L
  nn <- dims[2] - n1 # number of column sets
  nn1 <- nn - 1L
  nn2 <- 2^nn1
  # get all possible row values that will not lead to a transition back to the
  # initial state
  m <- as.matrix(expand.grid(rep(list(0:1), nn1)))
  m <- unique(
        function(i) {
          mm <- matrix(1L, nn2, dims[2])
          if (i > 1L) mm[, 1:(i - 1L)] <- m[, 1:(i - 1L)]
          if (i + n <= dims[2]) mm[, (i + n):dims[2]] <- m[, i:nn1]
  # for each row in m, which column sets would it increment?
  # aggregate unique counts to feed the creation of a transition matrix
  dtM <- as.data.table(
    m %*% (diag(0L, dims[2], nn) + rep(c(rep(1L, n), integer(nn)), nn)[1:(dims[2]*nn)]) == n
    , .(k = .N), eval(paste0("V", 1:nn))
  mTrans <- as.matrix(dtM[, 1:nn])
  class(mTrans) <- "numeric"
  # build a data.table listing all the state changes that don't end in the
  # initial state, and the number of ways the state change can occur
  dtTrans <- setorder(
        paste0("from", 1:nn)
      mTrans, dtM$k, n, nn
    , `:=`(
      state1 = rleidv(dtTrans[, 1:nn]),
      state2 = setorder(
          dtTrans[, (nn + 1L):(2L*nn)][, idx := .I]
          , state2 := rleidv(.SD) + 1L, .SDcols = 1:nn
        ], idx
  # build the transition matrix, calculate dims[1] transitions, and get the
  # probability of ending in the absorbing state (when the desired submatrix has
  # occurred)
  nStates <- max(dtTrans$state2)
  m <- matrix(0, nStates, nStates)
  m[as.matrix(dtTrans[, state1:state2])] <- dtTrans$k
  m[nStates, nStates] <- 2^dims[2] # absorbing state
  m[, 1] <- 2^dims[2] - sum(dtM$k)
  m[nStates, 1] <- 0
  matpow(m*p^dims[2], dims[1])[1, nStates]

# test the pGridInArray function
pGridInArray(2L, c(5L, 3L))
#> [1] 0.3262634
findSquare(as.matrix(expand.grid(rep(list(0:1), 15))), 5L, 3L, 2L)/2^15
#> [1] 0.3262634

pGridInArray(3L, c(5L, 4L))
#> [1] 0.009971619
findSquare(as.matrix(expand.grid(rep(list(0:1), 20))), 5L, 4L, 3L)/2^20
#> [1] 0.009971619

pGridInArray(5L, c(100L, 5L))
#> [1] 2.772544e-06
#> [1] 2.772544e-06

# 5-by-5 filled submatrix in 100 rows x 5-11 columns
p5 <- vapply(5:11, function(i) pGridInArray(5L, c(100L, i)), 0)
dt <- data.table(
  columns = 5:11,
  prob = p5,
# probability of finding the submatrix in the last 5 columns
# given that it is not anywhere else
  marg = c(NA, diff(p5)/(1 - first(p5, -1)))
# project the probability of finding the submatrix in 100 columns
# using the last "marginal" probability
  , proj := 1 - (1 - prob)*(1 - marg)^(100 - columns)
#>   columns         prob         marg         proj
#> 1       5 2.772544e-06           NA           NA
#> 2       6 5.454110e-06 2.681573e-06 0.0002574892
#> 3       7 8.135459e-06 2.681364e-06 0.0002574695
#> 4       8 1.081678e-05 2.681347e-06 0.0002574679
#> 5       9 1.349810e-05 2.681345e-06 0.0002574678
#> 6      10 1.617941e-05 2.681346e-06 0.0002574678
#> 7      11 1.886071e-05 2.681346e-06 0.0002574678

# show that the estimate is converging
diff(last(dt$proj, -1))
#> [1] -1.969634e-08 -1.576093e-09 -1.256927e-10  2.543288e-11 -5.694340e-13

# final estimate
#> [1] 0.0002574678

Unfortunately, I don't have a method to place error bounds on the final estimate, nor can I prove that the estimate is strictly converging with each iteration.


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.