Preliminary Survey :

N <- 100
M = 900
Mbar <- M/N

f <- c(3/6, 2/6, 1/2, 1/2, 2/2, 1/2,  1/2, 0/2, 6/8, 0/2, 0/2, 8/8, 7/8, 4/4, 4/4, 4/4, 2/4, 1/6, 1/10, 2/4, 6/10, 6/10, 4/8, 8/12)


mi <- c(6, 6, 2, 2, 2, 2, 2, 2, 8, 2, 2, 8, 8, 4, 4, 4, 4,6,10,4,10, 10,8,12)
ai <- c(3, 2, 1, 1, 2, 1, 1, 0, 6, 0, 0, 8, 7, 4,4,4,2,1,1,2,6,6,4,8) 
mbar <- mean(mi)

phat_per <- ai/mi
phat <- mean(phat_per)

# est variance  :
sp2 <- sum((ai - phat*mi)^2)/(length(mi)-1)
B <- .1 # off by 5 pct : bound-of-err
D <- ((B^2)*mbar^2)/4
n_est <- (N*sp2)/(N*D + sp2)
n_est <- ceiling(n_est) #round to nxt whole

Drawing Sample :

set.seed(123)
d1_samp <- sample(1:100, 30)

d1_samp |> as.data.frame() # True Order 
##    d1_samp
## 1       31
## 2       79
## 3       51
## 4       14
## 5       67
## 6       42
## 7       50
## 8       43
## 9       97
## 10      25
## 11      90
## 12      69
## 13      57
## 14       9
## 15      72
## 16      26
## 17       7
## 18      95
## 19      87
## 20      36
## 21      78
## 22      93
## 23      76
## 24      15
## 25      32
## 26      84
## 27      82
## 28      41
## 29      23
## 30      27
d1_samp |> as.data.frame() |> dplyr::arrange(d1_samp)
##    d1_samp
## 1        7
## 2        9
## 3       14
## 4       15
## 5       23
## 6       25
## 7       26
## 8       27
## 9       31
## 10      32
## 11      36
## 12      41
## 13      42
## 14      43
## 15      50
## 16      51
## 17      57
## 18      67
## 19      69
## 20      72
## 21      76
## 22      78
## 23      79
## 24      82
## 25      84
## 26      87
## 27      90
## 28      93
## 29      95
## 30      97