%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SCRIPT ANMT_Lec2_binomialLikelihoodTableWithPlots
%
% Simulates Binomial random samples from Binomial(.3), and determines the
% log-likelihood of each given a vector of guesses 0.01:0.01:.99.
%
% Instead of showing the log-likelihoods in a table, it displays 3 plots.
%
% INPUT: Random samples from Binomial(.3). Created within program.
%
% OUTPUT: likelihoodTable(Grows dinamically): Table with likelihood values
%
% Cesar Zamudio, 8.31.2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Clean house
clear; clc; format bank;
%Simulate required random variables in three vectors.
%Remember that the whole sample is not needed when dealing with the
%binomial log-likelihood: only the SUM of the successes is needed. This is
%what is reported by the binornd function and is labeled "x" in Assig. 2.
bin50=binornd(50,.3);
bin100=binornd(100,.3);
bin1000=binornd(1000,.3);
%Log-likelihood table
logLikelihoodTable=zeros(3,5);
%%%%%%%%%%%%%%%%%%%%%%%%%
% ALGORITHM %
%%%%%%%%%%%%%%%%%%%%%%%%%
%Obtain log-likelihoods for first sample (N=50)
N=50;
%The cycle's "i" index will be used to represent the parameter p.
%The only problem with using "i" as p is that we cannot use i
%directly to access our likelihood table matrix cells (because there's
%no such thing as a location in, say, (1, 1.5). So, let's create an
%index that will grow along with i, to know exactly where we are.
index=1;
%Finally, understand that "index" is started on the OUTSIDE of the cycle.
%This is so it will not be reset WITHIN the cycle, thus causing the index
%to never advance. On the other hand, p is started INSIDE of the
%cycle, because we want it to change as i changes.
for i=0.01:0.01:.99
p=i;
%For each value of p, we must calculate the log-likelihood. The
%log-likelihood has three parts, so we will calculate them separately,
%and then sum them to form the final log-likelihood.
%IMPORTANT: In general, the constant is irrelevant for the
%log-likelihood calculation. We will calculate the first part of the
%log-likelihood, but because of its size, it's best to omit it.
%PART 1. Obtain the log of constant k as defined in Assignment 2.
%This is a straightforward calculation.
LOGLIKPART1=factorial(N)/ ( factorial(bin50)*factorial(N-bin50) ) ;
%PART 2. x*log(p). Also straightforward.
LOGLIKPART2=bin50*log(p);
%PART 3. (n-x)*log(1-p). Also a simple calculation.
LOGLIKPART3=(N-bin50)*log(1-p);
%Finally, form the log-likelihood.
LOGLIKELIHOOD=LOGLIKPART2+LOGLIKPART3;
%Store the log-likelihood in the corresponding table cell.
logLikelihoodTable(1,index)=LOGLIKELIHOOD;
%MAKE SURE TO INCREASE YOUR INDEX, so you update the position where you
%are.
index=index+1;
end
%Let's plot the log-likelihoods for each value of p we used above.
%The "plot" function takes the values of X, of Y, and then plots them.
%In our case, X=p, the parameter we're assessing
% Y=The likelihood we calculated for each value of p.
%We will have a plot for each of our three samples.
subplot(1,3,1)
%REMEMBER TO HOLD ON; OTHERWISE YOU WILL OVERWRITE YOUR PLOT
hold on
%Let's plot the first sample results
plot((0.01:0.01:.99),logLikelihoodTable(1,:));
%Let's put some labels on the plot
title('Binomial Log-Likelihood plot, N=50, p=.3')
xlabel('Values of p')
ylabel('Log-likelihood score')
%The above algorithm is simply copied and updated to deal with the other
%samples. Thus, we will not comment this section. Realize that there is a
%much more efficient way to do this: have TWO cycles, one that creates
%random Binomial data, and another one that calculates the likelihoods.
%However both accomplish the same thing.
%Algorithm for N=100
%Obtain log-likelihoods for second sample (N=100)
N=100;
index=1;
for i=0.01:0.01:.99
p=i;
%Just omit Part 1
LOGLIKPART2=bin100*log(p);
LOGLIKPART3=(N-bin100)*log(1-p);
LOGLIKELIHOOD=LOGLIKPART2+LOGLIKPART3;
%REMEMBER TO STORE VALUES IN THE NEXT ROW OF YOUR LIKELIHOOD TABLE
logLikelihoodTable(2,index)=LOGLIKELIHOOD;
%MAKE SURE TO INCREASE YOUR INDEX, so you update the position where you
%are.
index=index+1;
end
%Access NEXT subplot panel
subplot(1,3,2)
%Let's plot the second sample results
%REMEMBER TO ACCESS THE NEXT ROW OF YOUR LIKELIHOOD TABLE
plot((0.01:0.01:.99),logLikelihoodTable(2,:));
%Let's put some labels on the plot
title('Binomial Log-Likelihood plot, N=100, p=.3')
xlabel('Values of p')
ylabel('Log-likelihood score')
%Algorithm for N=1000
%Obtain log-likelihoods for second sample (N=1000)
N=1000;
index=1;
for i=0.01:0.01:.99
p=i;
%Just omit Part 1
LOGLIKPART2=bin1000*log(p);
LOGLIKPART3=(N-bin1000)*log(1-p);
LOGLIKELIHOOD=LOGLIKPART2+LOGLIKPART3;
%REMEMBER TO STORE VALUES IN THE NEXT ROW OF YOUR LIKELIHOOD TABLE
logLikelihoodTable(3,index)=LOGLIKELIHOOD;
%MAKE SURE TO INCREASE YOUR INDEX, so you update the position where you
%are.
index=index+1;
end
%Access NEXT subplot panel
subplot(1,3,3)
%Let's plot the second sample results
%REMEMBER TO ACCESS THE NEXT ROW OF YOUR LIKELIHOOD TABLE
plot((0.01:0.01:.99),logLikelihoodTable(3,:));
%Let's put some labels on the plot
title('Binomial Log-Likelihood plot, N=500, p=.3')
xlabel('Values of p')
ylabel('Log-likelihood score')