%find_the_best
% Model of the classic 'secretary' problem. Described very well in
% Impossible by Julian Havil pp169-171.
%
% Idea is you have n candidates invited for an interview. You select k at
% random and determine who the best is within this subset. All are
% rejected. You then choose the first candidate from the remaining n-k who
% is better than all those candidates seen so far.
%
% What should k be to maximize the chance that this strategy yields the
% best candidate overall?
%
% The answer is k = round( (n-1)/e ) i.e. about 0.37% of n-1, which results
% in a maximum probability of choosing the best (as n becomes large) of
% 1/e.
%
% LAST UPDATED by Andy French May 2021

function find_the_best

%Number of candidates
n = [10,50,100,150]; 

%Graph fontsize
fsize = 18; 

%

%Plot Pbest vs k for each n
figure('name','find the best','color',[1 1 1]);
hold on; grid on; box on; set(gca,'fontsize',fsize);
xlabel('Rejected candidates k')
ylabel('P(best chosen after rejecting k)');
title({'Probability of best candidate chosen from','n candidates after rejecting the first k'});
rgb = rand(length(n),3);
for i=1:length(n)
    [Pe,k] = Pbest_exact( n(i) );
    plot(k,Pe,'linestyle','-','color',rgb(i,:));
    lgnd{i} = ['n = ',num2str(n(i))];
end
xlimits = get(gca,'xlim'); plot(xlimits,exp(-1)*[1,1],'k-');
for i=1:length(n)
    [Pa,k] = Pbest_approx( n(i) );
    [Pe,k] = Pbest_exact( n(i) );
    plot(k,Pe,'marker','.','color',rgb(i,:));
    plot(k,Pa,'linestyle','--','color',rgb(i,:));
    plot( (n(i)-1)/exp(1), exp(-1),'*','color',rgb(i,:) );
end
legend(lgnd); 

%Print PNG and close
print(gcf,'best candidate proability.png','-dpng','-r300');
close(gcf);

%%

%Harmonic sum
function y = H(n)
y = zeros(size(n));
for k=1:numel(n)
    nn = 1:n(k); y(k) = sum( 1./nn );
end

%%

%Exact solution for probability of best candidate chosen from n candidates
%after rejecting the first k
function [P,k] = Pbest_exact(n)
k = 0:n-1;
P = (1/n)*( 1 + k.*( H(n-1) - H(k) ) );

%%

%Approximate solution for probability of best candidate chosen from n candidates
%after rejecting the first k
function [P,k] = Pbest_approx(n)
k = 0:n-1;
P = (1/n)*( 1 + k.*log( (n-1)./k ) );

%End of code