% L4A demo: ML-estimation for proportion of defectives

n = 200;
k = 22;

% Try all these values of p
plist = linspace(0.01, 0.99, 1000)';

% this is just to allocate space for the results
L = zeros(size(plist));
l = zeros(size(plist));

% this is just to turn off Matlab's warning that
% the nchoosek results are not exact.  They are still
% "accurate to 15 digits" so good enough for us.
warning('off','MATLAB:nchoosek:LargeCoefficient');

for i=1:length(plist)
    p = plist(i);
    L(i) = nchoosek(n,k) * p^k * (1-p)^(n-k);
    l(i) = log(L(i));
end

[lmax,imax] = max(l);

subplot(2,1,1)
plot(plist,L,'r-')
hold on
plot([plist(imax) plist(imax)], [0 L(imax)], 'ko--', 'markerfacecolor', 'k')
grid on
ylabel('likelihood')
xlabel('p')

fprintf('Best p that we found (from grid search) was %.6f\n', plist(imax));

subplot(2,1,2)
plot(plist,l,'r-')
hold on
plot([plist(imax) plist(imax)], [-500 l(imax)], 'ko--', 'markerfacecolor', 'k')
grid on
ylabel('log likelihood')
xlabel('p')
set(gca, 'ylim', [min(l) max(50, max(l))]);