EDIT: based on others' feedback.
Here's an interesting thought:
set the list of {1,2,3,4,5,6,7}. Throw the coin for each element in the list sequentially. If it lands head side up for a particular element, remove the number from the list. If all the numbers from a particular iteration of the list are removed, repeat the sampling. Do so until only one number remains.
drop.one <- function(x, p) {
drop <- runif(length(x)) < p
if (all(drop))
return(x)
return(x[!drop])
}
sample.recur <- function(x, p) {
if (length(x) > 1)
return(sample.recur(drop.one(x, p), p))
return(x)
}
# x <- c(1:7,7:1)
x <- 1:7
p <- 0.01
out <- replicate(1e5, sample.recur(x, p))
round(prop.table(table(out)), 2)
gives me an approximately uniform distribution
> round(prop.table(table(out)), 2)
out
1 2 3 4 5 6 7
0.14 0.14 0.15 0.14 0.14 0.14 0.14
It's interesting to note (if I haven't made a dire mistake) that this produces a different result than generating $N$ binomial outcomes as the sum of 13 tosses of the coin (counting 0 heads as an outcome) and mapping the {0,1,2,...,12,13} index onto the earlier list of {1,2,3,...,3,2,1}. I don't quite know how to prove that my method works.
##Evaluation of expectation value for number of coin throws
The expectation value for the number coin throws can be calculated using the transition matrix below (answering the question when we start with $x$ non-eliminated numbers then what is the probability to get to $y$ non-eliminated numbers)
$$M = \begin{bmatrix} q^7 & 0 & 0 & 0 &0 & 0 & 1 & 1\\ 7p^1q^6 & q^6 & 0 & 0 & 0 & 0 & 0 & 0 \\ 21p^2q^5 & 6p^1q^5 & q^5 & 0 & 0 & 0 & 0 & 0\\ 35 p^3q^4 & 15 p^2q^4 & 5q^4 & q^4 & 0 & 0 & 0 & 0\\ 35 p^4q^3 & 20 p^3q^3 & 10 p^2q^3 & 4 p^1q^3 & q^3 & 0 & 0 & 0\\ 21p^5q^2 & 15 p^4q^2 & 10 p^3q^2 & 6 p^2q^2 & 3 p^1q^2 & q^2 & 0 & 0\\ 7p^6q^1 & 6 p^5q^1 & 5 p^4q^1 & 4 p^3q^1 & 3p^2q^1 & 2p^1q^1 & 0 & 0\\ p^7 & p^6 & p^5 & p^4 & p^3 & p^2 & 0 & 0 \end{bmatrix}$$
The eigenvector associated with the eigenvalue 1 (which can be found by solving $(M-I)v=0$) depicts how much time is relatively spend in what state. Then 7th state is how often you will be able to draw a number from 1 to 7. The other states tell how many coin flips it costs.
Below is the image which compares with the answer from NcAdams which has expectation value for coin throws being $E(n) = \frac{24}{7}p(1-p)$
Remarkable is that the method performs better roughly for $p>2/3$. But also the performance is non-symmetric. A symmetric and better overall performance could be made when a probabilistic switching rule would be made which changes the decision rule from tails to heads when heads happens to be improbable.
###Solution found with wxMaxima
M: matrix(
[(1-p)^7, 0, 0,0,0,0,1,1],
[7* p*(1-p)^6, (1-p)^6, 0,0,0,0,0,0],
[21*p^2*(1-p)^5, 6*p*(1-p)^5, (1-p)^5,0,0,0,0,0],
[35*p^3*(1-p)^4, 15*p^2*(1-p)^4, 5*p*(1-p)^4,(1-p)^4,0,0,0,0],
[35*p^4*(1-p)^3, 20*p^3*(1-p)^3, 10*p^2*(1-p)^3,4*p*(1-p)^3,(1-p)^3,0,0,0],
[21*p^5*(1-p)^2, 15*p^4*(1-p)^2, 10*p^3*(1-p)^2,6*p^2*(1-p)^2,3*p*(1-p)^2,(1-p)^2,0,0],
[7* p^6*(1-p)^1, 6*p^5*(1-p), 5*p^4*(1-p),4*p^3*(1-p),3*p^2*(1-p),2*(1-p)*p,0,0],
[p^7, p^6, p^5,p^4,p^3,p^2,0,0]
);
z: nullspace(M-diagmatrix(8,1));
x : apply (addcol, args (z));
t : [7,6,5,4,3,2,0,0];
plot2d(t.x/x[7],[p,0,1],logy);
Calculations in R
# plotting empty canvas
plot(-100,-100,
xlab="p",
ylab="E(total flips)",
ylim=c(10,1000),xlim=c(0,1),log="y")
# plotting simulation
for (p in seq(0.1,0.9,0.05)) {
n <- rep(0,10000)
for (i in 1:10000) {
success = 0
tests = c(1,1,1,1,1,1,1) # start with seven numbers in the set
count = 0
while(success==0) {
for (j in 1:7) {
if (tests[j]==1) {
count = count + 1
if (rbinom(1,1,p) == 1) {
tests[j] <- 0 # elliminate number when we draw heads
}
}
}
if (sum(tests)==1) {
n[i] = count
success = 1 # end when 1 is left over
}
if (sum(tests)==0) {
tests = c(1,1,1,1,1,1,1) # restart when 0 are left over
}
}
}
points(p,mean(n))
}
# plotting formula
p <- seq(0.001,0.999,0.001)
tot <- (7*(p^20-20*p^19+189*p^18-1121*p^17+4674*p^16-14536*p^15+34900*p^14-66014*p^13+99426*p^12-119573*p^11+114257*p^10-85514*p^9+48750*p^8-20100*p^7+5400*p^6-720*p^5)+6*
(-7*p^21+140*p^20-1323*p^19+7847*p^18-32718*p^17+101752*p^16-244307*p^15+462196*p^14-696612*p^13+839468*p^12-806260*p^11+610617*p^10-357343*p^9+156100*p^8-47950*p^7+9240*p^6-840*p^5)+5*
(21*p^22-420*p^21+3969*p^20-23541*p^19+98154*p^18-305277*p^17+733257*p^16-1389066*p^15+2100987*p^14-2552529*p^13+2493624*p^12-1952475*p^11+1215900*p^10-594216*p^9+222600*p^8-61068*p^7+11088*p^6-1008*p^5)+4*(-
35*p^23+700*p^22-6615*p^21+39235*p^20-163625*p^19+509425*p^18-1227345*p^17+2341955*p^16-3595725*p^15+4493195*p^14-4609675*p^13+3907820*p^12-2745610*p^11+1592640*p^10-750855*p^9+278250*p^8-76335*p^7+13860*p^6-
1260*p^5)+3*(35*p^24-700*p^23+6615*p^22-39270*p^21+164325*p^20-515935*p^19+1264725*p^18-2490320*p^17+4027555*p^16-5447470*p^15+6245645*p^14-6113275*p^13+5102720*p^12-3597370*p^11+2105880*p^10-999180*p^9+371000
*p^8-101780*p^7+18480*p^6-1680*p^5)+2*(-21*p^25+420*p^24-3990*p^23+24024*p^22-103362*p^21+340221*p^20-896679*p^19+1954827*p^18-3604755*p^17+5695179*p^16-7742301*p^15+9038379*p^14-9009357*p^13+7608720*p^12-
5390385*p^11+3158820*p^10-1498770*p^9+556500*p^8-152670*p^7+27720*p^6-2520*p^5))/(7*p^27-147*p^26+1505*p^25-10073*p^24+49777*p^23-193781*p^22+616532*p^21-1636082*p^20+3660762*p^19-6946380*p^18+11213888*p^17-
15426950*p^16+18087244*p^15-18037012*p^14+15224160*p^13-10781610*p^12+6317640*p^11-2997540*p^10+1113000*p^9-305340*p^8+55440*p^7-5040*p^6)
lines(p,tot)
#plotting comparison with alternative method
lines(p,3*8/7/(p*(1-p)),lty=2)
legend(0.2,500,
c("simulation","calculation","comparison"),
lty=c(0,1,2),pch=c(1,NA,NA))
