# Given a number of four-person families, how many people need to be selected so 50% of the time, half of them are in familes with another selected?

Doing some worldbuilding and was thinking about how connected my world is. I had thoughts on dealing with how I could calculate pedigree collapse and such, but wanted to use a simpler model to think about for now, so:

Given a number of four-person families, how many people need to be selected so 50% of the time, half of them are in familes with another selected?

I know that it's obviously at least double the number of families, a.k.a half of the number of people, and less than everyone except two.

For 1 family, that's the trivial case. Any two people will do, and that's the minimum number for the selected to be in families with each other.

Suppose people are selected independently with probability $$p$$. The number of people $$M$$ selected in a family of four has a $$M\sim Binomial(4,p)$$ distribution, and you want $$P(M>1)/P(M>0)=1/2$$

In code

curve(pbinom(1,4,x,lower.tail=FALSE)/pbinom(0,4,x,lower.tail=FALSE),
ylim=c(0,1),ylab="P(another selected in family)",
xlab="proportion selected")
abline(h=1/2,lty=2,col="purple")


If 1/3 of the people are selected, your probability is just over 1/2. I might suspect there's a more elegant reason why it's 1/3, except that it isn't exactly 1/3. Numerical optimisation gives 0.32809 • I'd like a formula that outputs an integer number, but this looks useful, at least! Thanks! ... I guess if it's not an integer, you round up? 0.32809*(Number of People). If Number = 4, it should result in 2, so rounding up. Jun 8 at 3:35
• The ratio of probabilities is $$\frac{1-(1-p)^4-4p(1-p)^3}{1-(1-p)^4}$$ which leads to $p=(20-c-22/c)/21$, where $c=(694-42\sqrt{267})^{1/3}$, so indeed $p\simeq 0.32809<1/3$. Jun 8 at 4:45
• @Malady, there is no rounding — the binomial distribution describes the probabilities of 0,1,2,3, or 4 members of a family being chosen. Jun 8 at 4:47
• I believe your answer finds p such that P(# of families with 2 or more > # of families with 1) = 0.5. But I don't think that's what the question asks. I think it asks that P(# of individuals in families with 2 or more > # of individuals in families of 1) = 0.5 Jun 8 at 15:13

Here is a simulation of the sampling and calculation, done in R. If you don't have R installed you can copy the code below and run it here, https://rdrr.io/snippets/. I include two different ways of sampling.

1. Out of the $$N = num\_fams*4$$ individuals, define $$n = roundup(p*N)$$ and sample $$n$$ individuals (without replacement). This guarantees that exactly $$n$$ are sampled.
2. For each person, sample them with probability $$p$$. The mean number of people sampled is $$n = p*N$$, and there is some variation in the exact number sampled. The number sampled follows a binomial distribution.

My approach is to simulate a sample of $$N$$ people as 1s and 0s, where 1s indicate that person has been sampled. Add them up in groups of 4, giving num_fams values ranging from 0-4, then count the number of values that equal to 1, saved as num_singletons. That gives the number of families that had exactly 1 person sampled, and also equals the number of people who were not in a family with another selected.

Then I classify the simulation as a "success" in two ways. 1) if the number of singletons is less than half the total number sampled we count that as a success. 2) is inspired by @ThomasLumley's answer, if the number of singletons is less than the number of families with 2 or more sampled we count that as a success.

num_fams <- 200 # change as desired
num_sims <- 10000 # change as desired
constant <- (694 - 42*sqrt(267))^(1/3) # From CV answer
proportion <- (20 - constant - (22/constant))/21 # From CV answer

# fixed number to sample
num_to_sample <- ceiling(num_fams*4*proportion)
set.seed(1234)
successes <- rep(0, num_sims)
successful_fams <- rep(0, num_sims)
simulated_singletons <- rep(0, num_sims)
for(i in 1:num_sims) {
population <- rep(0, num_fams*4)
population[sample(1:length(population), num_to_sample)] <- 1
num_in_fams <- do.call(cbind,
lapply(seq(1, length(population), 4),
function(i) sum(population[i:(i+3)])))
num_singletons <- sum(num_in_fams == 1)
simulated_singletons[i] <- num_singletons
if((num_singletons/sum(population)) < 0.5) {
successes[i] <- 1
}
if(sum(num_in_fams > 1) > num_singletons) {
successful_fams[i] <- 1
}
}
mean(simulated_singletons)
#>  79.7053
mean(successes)
#>  1
mean(successful_fams)
#>  0.4815

# random number to sample
set.seed(1234)
successes <- rep(0, num_sims)
successful_fams <- rep(0, num_sims)
simulated_singletons <- rep(0, num_sims)
for(i in 1:num_sims) {
population <- rbinom(num_fams*4, 1, proportion)
num_in_fams <- do.call(cbind,
lapply(seq(1, length(population), 4),
function(i) sum(population[i:(i+3)])))
num_singletons <- sum(num_in_fams == 1)
simulated_singletons[i] <- num_singletons
if((num_singletons/sum(population)) < 0.5) {
successes[i] <- 1
}
if(sum(num_in_fams > 1) > num_singletons) {
successful_fams[i] <- 1
}
}
mean(simulated_singletons)
#>  79.6984
mean(successes)
#>  1
mean(successful_fams)
#>  0.4792


Created on 2022-06-08 by the reprex package (v2.0.1)

As we can see, taking $$p = 0.3280935$$ results in a roughly 50% chance of the second type of success, but is far too high for the first type of success. There are on average 79.7 individuals without any sampled family members, but we are sampling 262.4 individuals on average. So only ~30% of individuals are alone in their family.

Fiddling with the proportion, it looks like taking $$p = 0.208$$ makes the first type of success have a probability of about 50%.