1000人の女性の中の1位とパートナーになるには？その２．解析的に

某氏より、これは有名な秘書問題だよと教えていただいた。
https://ja.wikipedia.org/wiki/%E7%A7%98%E6%9B%B8%E5%95%8F%E9%A1%8C

これで簡単に解析的に解けると。

n=1000

f <- function(r){
      sum(1 / (r:n - 1)) * (r - 1) / n  
}
f2 <- Vectorize(f)
ans <- f2(2:1000)

which.max(ans)
max(ans)

簡単に最大となる368が出た。
36.81956　%となる

nmax <- 1000
tbl_h <- matrix(0, ncol=nmax, nrow=3)
rownames(tbl_h) <- c("N", "best i", "max p")

tbl_h[,1] <- c(1,0,0.5)
for(i in 2:nmax){
      f_temp <-  Vectorize(function(r){(sum(1/(r:i-1)) * (r - 1) / i)})  
      temp_1 <- f_temp(2:i)
      tbl_h[1, i] <- i
      tbl_h[2, i] <- which.max(temp_1)
      tbl_h[3, i] <- max(temp_1)
      
}

plot(tbl_h[2,1:10], type="p",xlab="N人の場合", ylab="i番目までを基準にすると最大確率") #はじめの10人
plot(tbl_h[3,], type="l",xlab="N人の場合", ylab="1位の女性とパートナーになれる最大確率")
abline(h = 1 / exp(1), col = "red" )

Nが増えると、1/eに近づいていくのが確認できた。