socks

socks.r

---
title: The mystery of Nick's unpaired socks
Authors: The Balfour Road consortium
---

## Introduction

Nick has only 16 socks and none of them form a pair. There are several mechanisms that could give rise to his predicament. The most likely follows a scenario approximating the following.

Nick buys a pair of socks, and after some period of time, he loses one of the socks. Nick then buys another pair of socks, and eventually loses one of those too. And so on. If the time to loss of a single sock is Poissonly distributed with $\lambda$ time to event, then the expected amount of time taken to reach a situation with $n$ unpaired socks is simply $n\lambda$. 

For example, if it takes Nick 3 months to lose a single sock, and he has 16 unpaired socks, he must be around 4 years old.

There are other scenarios that could give rise to Nick's predicament, and they are considered below.

## 1. Nick starts with 32 pairs of socks and loses one sock at a time

### 1a. Nick never wears unpaired socks

This scenario is much like that introduced at the beginning. Nick will only ever be able to lose socks that are paired because those are the only ones he wears. The scenario of being left with 16 unpairable socks is simply a matter of time. In order to get the correct estimated time we assume, as we did previously, that Nick only ever buys unique pairs of socks.

### 1b. Nick does wear unpaired socks

After having lost one sock, the next lost sock must not be the remaining sock of the pair from the first lost sock. After losing the second sock, the next lost sock must not be from either of the first two pairs lost, and so on. This can be expressed simply as follows:

$$
\begin{aligned}
p(no pairs) & = \prod^{n}_{i=2} \frac{2n-2i}{2n-i}
& = \frac{2n!!}{2n!}
\end{aligned}
$$

The probability of starting with 16 pairs of socks and being left with 16 unpairable socks is simply

```{r}

dfactorial <- function(x)
{
	start <- ifelse(x %% 2 == 0, 2, 1)
	return(prod(seq(start, x, 2)))
}

dfactorial(5) / factorial(5)



sim1 <- function(npair_start, n_left)
{
	socks <- rep(1:npair_start, each=2)
	nlost <- npair_start * 2 - n_left
	lost <- array(NA, nlost)
	for(i in 1:nlost)
	{
		ind <- sample(1:length(socks), 1)
		lost[i] <- socks[ind]
		socks <- socks[-ind]
	}
	return(any(duplicated(socks)))
}

sim1(16, 16)

sim <- function(npair_start, n_left, nsim)
{
	out <- array(0, nsim)
	for(i in 1:nsim)
	{
		out[i] <- sim1(npair_start, n_left)
	}
	return(1 - sum(out) / nsim)
}

poss <- 16:32
res <- rep(0, length(poss))
for(i in 1:length(poss))
{
	message(i)
	res[i] <- sim(poss[i], 16, 50000)
}

dat <- data.frame(res=res, poss=as.factor(poss))
library(ggplot2)
ggplot(dat, aes(x=poss, y=res)) +
geom_bar(stat="identity") +
geom_label(aes(y=res+0.005, label=round(res, 3))) +
labs(x="Starting number of pairs", y="Probability of 16 unmatching socks")


theor <- function(npair_start, nlost)
{
	nlost <- min(npair_start / 2, nlost)

	if(nlost < (npair_start *2)) return(1)

	l <- array(1, nlost)
	for(i in 2:(nlost))
	{
		l[i] <- ((2*npair_start - 2*i) / (2*npair_start - i))
	}
	print(l)
	return(prod(l))
}

theor(16, 32)

```