set.seed(2015)

simulate_seats <- function(seats, simulations) {
  m <- matrix(seq_len(seats), nrow = seats, ncol = simulations)

  for (i in seq(seats - 1)) {
    if (i > 1) {
      # take someone's seat if this is not the right one
      taken <- which(m[i, ] != i)
    } else {
      # replace all seats
      taken <- seq_len(simulations)
    }

    # Each remaining person switches seats with someone remaining
    switch_with <- sample(seq(i, seats), length(taken), replace = TRUE)

    replacements <- m[cbind(switch_with, taken)]
    m[cbind(switch_with, taken)] <- m[i, taken]
    m[i, taken] <- replacements
  }
  m
}

# perform 1,000,000 simulations with 100 seats
sim <- simulate_seats(100, 1000000)

The probability of the final (100th) person sitting in the correct seat is:

mean(sim[100, ] == 100)
## [1] 0.498978

That is, 49.9%.

We can also look at the % of sitting in the correct seat as a function of the person’s assigned seat.

correct_by_seat <- rowMeans(sim == seq_len(100))

plot(correct_by_seat, type = "l",
     xlab = "Assigned seat",
     ylab = "Probability of correct seat")

Here we see that the first person has (of course) a 1/100 chance of sitting in the right seat, but that after that most people are likely to sit in the correct one, until the last few to get on the plane.