Repeated catch–mark–release (urn problem) Imagine I have a small pond with some fish in it and I want to estimate its population size. I don't have a lot of resources at hand, in fact all I do have is a fishing pole, pen and paper, and a lot of patience. 
One morning I sit down at the end of the pier and start fishing. For every fish I pull I examine it for a mark. If it has none I write down a 0 in my note book, give the fish a mark, and release it. If it already has a mark I write down a 1 in my note book and release the fish. I repeat this procedure 40 times and at the end have the following sequence: 
0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1 0 0 1 0 0 0 1 0 0 0 1 0 1

How can I use this data to estimate the population size of the fish in my pond? And compared to a standard mark–recapture procedure where I capture just as many fish, in what ways are my method worse?
I assume ideal conditions, like constant population size and each fish being equally likely to get caught. 
So far all I have are the reasonably obvious. The probability of capturing a marked fish increase every time I capture a non-marked fish, which sounds a bit like a poisson process (no pun intended).
My problem seems also to be similar to this question, but different enough for me not to know what to do about it.
For any fellow R fiends and others interested, the above sequence was generated like this:
set.seed(1)
pool <- rep(0, 64)
capture <- vector()

for (i in 1:40) {
capture[i] <- sample(pool, 1)
if (capture[i]==0) pool <- c(pool[-1], 1)
}

 A: Let $x_i=1$ denote the event that the $i$'th fish caught is already marked and $x_i=0$ the event that you catch a new fish which you then mark.  Prior to catch $i$, the number of fish marked (the cumulative number of zeros in your sequence) is 
$$
\sum_{j=1}^{i-1} (1-x_j)=i-1-\sum_{j=1}^{i-1}x_j.
$$
Hence, letting $N$ denote the unknown population size, the probability mass function of $x_i$ conditional on past observation is
$$
p(x_i|x_1,\dots,x_{i-1}) = f(x_i,\frac{i-1-\sum_{j=1}^{i-1}x_j}N)
$$
where $f(x,p)$ is the probability mass function of a Bernoulli variable with parameter $p$ (dbinom(  ) in the code below).  By the probability chain rule, the total likelihood is
$$
L(N) = p(x_1,\dots,x_n)=p(x_1)p(x_2|x_1)\cdots p(x_n|x_1,\dots,x_{n-1}).
$$
The following R code implements this and produces the plot below.
loglikelihood <- function(N,x) {
  n <- length(x) 
  i <- 1:n
  marked <- i - 1 - c(0,cumsum(x))[-n]
  sum(dbinom(x, size=1, prob=marked/N, log=TRUE))
}
NN <- 1:300
ll <- vector()
for (i in 1:length(NN)) 
  ll[i] <- loglikelihood(NN[i],x=capture)
plot(NN,exp(ll),type="l",ylab="Likelihood",xlab="Population size N")


