%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SCRIPT Main_LL_Normal
% 
% This function runs MLE Estimation for a simulated Normal distribution.
% Estimates mu and sigma.
%
% INPUT: DATA - Simulated Normal data.
%
% OUTPUT: MLE estimates of the Normal. Should recover them.
%
%                                                   César Zamudio 9.16.2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc; format bank;

%True parameters
mu=35;
sigma=7.2;

%Sample size
N=100000;

%Draw a sample of N normals
DATA=mu+sigma*randn(N,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAXIMUM LIKELIHOOD ESTIMATION ROUTINE %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Prepare minimizer options. Investigate these in MATLAB documentation.
options= optimset('Display', 'iter','TolX', 1e-16); 

%STARTING VALUES
%---------------
%This vector ("parameters") holds the starting values for the parameters.
%You are interested in two parameters (mu and sigma). Therefore, you need
%two parameters, which we'll initialize at 1.
STARTING_VALUES=ones(2,1);

%MINIMIZATION ROUTINE
%--------------------
%We will use function "fminunc". It is a minimization algorithm. USE THE
%SAME SYNTAX AS HERE in your own work. The pieces you're gonna work with
%are:
%STARTING_VALUES: The starting values you defined above.
%PARAMETERS: This variable is not defined yet. It will hold your estimated
%parameters (MLE estimates).
%LOG_LIKELIHOOD_NORMAL: Here you will put the code for whatever function
%you want to minimize. In this case it's the normal. The parameters it
%should receive are:
%1) x: These are the PARAMETERS - they only have a different name.
%2) DATA: That's your data.
%3) (...) You can pass on any number of parameters as soon as you defined
%that in your LOG_LIKELIHOOD_NORMAL function (or any function)

[PARAMETERS,fval,exitflag,output,grad,hessian2]=fminunc(@(x) LOG_LIKELIHOOD_NORMAL(x,DATA),STARTING_VALUES,options);

disp('Estimate of mu')
PARAMETERS(1,1) %Display MLE estimate of mean 
disp('Estimate of sigma2')
PARAMETERS(2,1) %Display MLE estimate of variance
disp('We are interested in sigma, not sigma square, so take square root:')
sqrt(PARAMETERS(2,1))

%Matlab's function MLE applied to the same data should yield the EXACT same
%results.
disp('MATLABs MLE and our MLE should be identical! Minor differences')
disp('can be observed due to rounding or minimization settings')
mle(DATA)
