Generate a matrix similar to what Brandon is describing:

bin_seq <- function(n, p) {
  x <- sample(c('1', '0'), size = n, prob = c(p, 1-p), replace = TRUE)
  paste(x, collapse = '')
}

P <- 0.3
N <- 8
x <- sapply(1:(N*N), function(i) bin_seq(N,P))
m <- matrix(x, nrow = N)

Okay, disregard the matrix. It provides no information. We just care about the proportion of ones and zeroes and the number of trials.

X <- paste(x, collapse = '')
num_trials <- nchar(X)
successes <- sum(as.numeric(str_split(X, '')[[1]]))

Now a Bayesian grid search. First create a grid.

n_grid <- 1000
g <- seq(0, 1, length.out = n_grid)

Now compute the log-likelihood.

ll <- dbinom(successes, size = num_trials, prob = g, log = TRUE)

Now using a flat prior. Add them together to create an unstandardized posterior. Then standardized it.

un_post <- ll + rep(1/n_grid, n_grid)
un_post <- exp(un_post)
post <- un_post / sum(un_post)

Now plot it.

qplot(g, post, geom = 'line', xlab = 'p', ylab = 'Posterior probability')

Sample from the distribution and create a credible interval.

library(tidybayes)
s <- sample(g, 1000, replace = T, prob = post)
ci <- hdi(s)
print(ci)
##           [,1]      [,2]
## [1,] 0.2702703 0.3453453