Question: The Secretary Problem

I was wondering if anyone could explain the theoretical reason (ie derive the
expression for the probability) why the optimal solution for the secretary problem
with n=100 is equal to 100/exp(1) and the probability is (1/exp(1))*100

Robert Israel tend to be good at these kind of things....:-)

The problem can be stated as follows (wikipedia):

 1. There is a single secretarial position to fill.
 2. There are n applicants for the position, and the value of n is known.
 3. The applicants can be ranked from best to worst with no ties.
 4. The applicants are interviewed sequentially in a random order, with each order being equally likely.
 5. After each interview, the applicant is accepted or rejected.
 6. The decision to accept or reject an applicant can be based only on the relative ranks of the applicants interviewed so far.
 7. Rejected applicants cannot be recalled.
 8. The objective is to select the best applicant. The payoff is 1 for the best applicant and zero otherwise.

 
The objective in the problem is to select the single best applicant. The optimal policy is to interview and reject the first n/e applicants (where e is the base of the natural logarithm) and then to accept the next who is better than the previous interviewed candidates. As n gets larger, the probability of selecting the best applicant from the pool goes to 1/e, which is around 37%. Whether one is searching through 100 or 100,000,000 applicants, the optimal policy will select the single best one about 37% of the time.


This can be seen by running the below code (it is a bit slow due to the large nsim number):


restart:
with(plots):
with(Statistics):
randomize():
with(LinearAlgebra):

X := proc (w) local nsim, n, s, j, i, sel, data;

nsim := 2000:
n := 100:

s := Transpose(Matrix([seq(Shuffle([seq(i, i = 1 .. n)]), i = 1 .. nsim)])):

for j to nsim do
for i to n do

if w < i and max(s[1 .. i-1, j]) < s[i, j] then sel[j] := s[i, j]; break else sel[j] := 0 end if

end do:
end do:


for j to nsim do

if sel[j] = max(s[1 .. n, j]) then data[j] := 1 else data[j] := 0 end if

end do:

evalf(100*add(data[j], j = 1 .. nsim)/nsim, 4) ;

end proc:

evalf(100/exp(1), 4);
100*evalf(1/exp(1), 4);

display({plot([seq([w, X(w)], w = 0 .. 100)]), plot([36.79, y, y = 0 .. 50], color = black, thickness = 3, linestyle = dash), plot(37, x = 0 .. 50, color = black, thickness = 3, linestyle = dash)});