---
title: "Mark and recapture likelihood"
author: ""
date: ""
output: html_document
---

Here is an illustration of how well the postulated binomial distributions
approximate the the actual relative frequency of counts of the number of
recaptured marked animals ($R$) and the number of unmarked captures ($Z$).

## Simulating the experiment

```{r}
set.seed(42)
N <- 1000
M <- 100
theta <- 0.05
```

Assume the true unknown population size is N=`r N` and we have marked M=`r M` at
a previous visit. Assume the catch probability is $\theta=`r 100*theta`$% 
(important that it is low i.e. $\theta$<<1). I.e. when we return to the population we
catch a random number of animals $K$ which must be approximately binomial with
parameters $N$ and $\theta$ (ignoring that we catch without replacement).

A single experiment could be performed in R like this:
```{r}
pop <- factor(c(rep(0, N-M), rep(1,M))) # Population: 0 unmarked, 1 marked
K <- rbinom(1, size = N, prob = theta) # Random number of catches
s <- sample(pop, size = K) # Single sample from population
tab <- table(s) # Table of counts of "0" and "1"
tab
```

We can then repeat this experiment many (`nrep`) times. In R the `replicate()`
function can be used to run the same code many times:

```{r}
nrep <- 10000
tab <- replicate(nrep, {
  K <- rbinom(1, size = N, prob = theta) # Random number of catches
  s <- sample(pop, size = K) # Single sample from population
  table(s) # Table of counts of "0" and "1"
})
```

The first 10 results for unmarked/marked:
```{r}
tab[,1:10]
```

Relative frequency of the number of recaptures compared with a binomial density
with parameters $M$ and $\theta$ (observed relative frequency as a black bar,
approximate binomial probability as red circle):

```{r}
Rtable <- table(Z=tab[2,])
R <- as.numeric(names(Rtable))
Rfreq <- as.numeric(Rtable)
plot(R, Rfreq/nrep, type = "h")
points(R, dbinom(R, size = M, prob = theta), col = 2)
```

Relative frequency of the number of unmarked captures compared with a binomial
density with parameters $U=N-M$ and $\theta$:

```{r}
Ztable <- table(Z=tab[1,]) # Count
Z <- as.numeric(names(Ztable))
Zfreq <- as.numeric(Ztable)
plot(Z, Zfreq/nrep, type = "h")
points(Z, dbinom(Z, size = N-M, prob = theta), col = 2)
```