Cool from an educational point of view, but also severely overkill: you can get the posterior distribution of how many socks someone has, given that they just ended up with 11 out of 11 unmatched socks, with SciPy and a little loop:
posterior = [hypergeom.pmf(11, M=n_pairs*2, n=n_pairs, N=11) * poisson.pmf(n_pairs, mu=20) for n_pairs in range(50)]
You don't need Metropolis or Approximate Bayesian Computation because there's maybe only 50 plausible answers, meaning you can just run through them all.
