%binomialExample.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reproduce Binomial PDFs in Figure 1 (Myung, 2003; p. 91) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1), clf;
p02 = binopdf(0:10, 10, 0.2); %PMF for all possible Ncorrect given N=10, p=0.2
subplot(2,1,1), bar(0:10,p02);
xlabel('Data y'); ylabel('f(y|n=10, w=0.2)');
text(9,0.3,'w = 0.2');
p07 = binopdf(0:10, 10, 0.7); %PMF for all possible Ncorrect given N=10, p=0.7
subplot(2,1,2), bar(0:10,p07);
xlabel('Data y'); ylabel('f(y|n=10, w=0.7)');
text(9,0.3,'w = 0.7');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reproduced Likelihood function in Figure 2 (Myung, 2003; p. 92) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 10; %number of Bernoulli trials in experiment
y = 7;  %single data observation (Ncorrect);
w = [0:0.01:1]; %range of parameter values (independent Bernoulli probability)
L = binopdf(y,N,w); %compute Likelihood function for Binomial model
figure(2), clf, plot(w,L);
xlabel('parameter w');
ylabel('Likelihood L(w|n=10, y=7)');
title('Binomial Likelihood Example');
%
%simple Maximum Likelihood Estimation for example above
%first, search 'L' for location of the maximum likelihood value
maxL = max(L);
%then find the value of 'w' associated with the maximum L
bestfit_w = w(find( L == maxL));
str1=sprintf('max(L) given by w = %04.2f', bestfit_w);
text(0.1,0.3,str1);

%some miscellaneous explorations of Myung(2003)
%log-likelihood approach (page 93)
N = 10; %number of Bernoulli trials in experiment
y = 7;  %single data observation (Ncorrect);
w = [0:0.01:1]; %range of parameter values (independent Bernoulli probability)
%compute log-likelihood (computational advantages of sum & products w/o
%requirements for exponentiation, etc.
%i.e., lnL = ln(N!/(y!*(N-y)!)) + y*ln(w) + (N-y)*ln(1-w)
lnL = log(factorial(N)/(factorial(y)*factorial(N-y)))+(y*log(w))+(3*log(1-w));
%second derivative
%d_lnL = y - (N*w)/(w*(1-w));
