%% Assess Model Specifications Using the Likelihood Ratio Test
% Compare two model specifications for simulated education and income data. The unrestricted model has the following loglikelihood:
% l(\beta,\rho) = -n\log\Gamma(\rho) + \rho\sum^n_{k=1}\log\beta_k + (\rho-1)\sum^n_{k=1}\log(y_k)-\sum^n_{k=1}y_k\beta_k,
% where
% \beta_k=\frac{1}{\beta+x_k}.
% x_k is the number of grades that person k completed.
% y_k is the income (in thousands of USD) of person k.
% That is, the income of person k given the number of grades that person k completed is Gamma distributed with shape \rho and rate \beta_k. The restricted model sets \rho = 1, which implies that the income of person k given the number of grades person k completed is exponentially distributed with mean \beta + x_k.
% The restricted model is H_0: \rho = 1. Comparing this model to the unrestricted model using lratiotest requires the following:
% The loglikelihood function
% The maximum likelihood estimate (MLE) under the unrestricted model
% The MLE under the restricted model
% Load the data.

load Data_Income1
x = DataTable.EDU;
y = DataTable.INC;

%
% To estimate the unrestricted model parameters, maximize l(\rho,\beta) with respect to \rho and \beta. The gradient of l(\rho,\beta) is
% \frac{\partial l(\rho,\beta)}{\partial \rho} = -n\psi(\rho) + \sum^n_{k=1}\log(y_k\beta_k)
% \frac{\partial l(\rho,\beta)}{\partial \beta} = \sum^n_{k=1}\beta_k(\beta_k y_k-\rho),
% where \psi(\rho) is the digamma function.

nLogLGradFun = @(theta) deal(-sum(-gammaln(theta(1)) - ...
    theta(1)*log(theta(2) + x) + (theta(1)-1)*log(y) - ...
    y./(theta(2)+x)),...
    -[sum(-psi(theta(1))+log(y./(theta(2)+x)));...
    sum(1./(theta(2)+x).*(y./(theta(2)+x)-theta(1)))]);

% nLogLGradFun is an anonymous function that returns the negative loglikelihood and the gradient given the input theta, which holds the parameters \rho and \beta, respectively.
%
% Numerically optimize the negative loglikelihood function using fmincon, which minimizes an objective function subject to constraints.

theta0 = randn(2,1); % Initial value for optimization
uLB = [0 -min(x)];   % Unrestricted model lower bound
uUB = [Inf Inf];     % Unrestricted model upper bound
options = optimoptions('fmincon','Algorithm','interior-point',...
    'FunctionTolerance',1e-10,'Display','off',...
    'SpecifyObjectiveGradient',true); % Optimization options

[uMLE,uLogL] = fmincon(nLogLGradFun,theta0,[],[],[],[],uLB,uUB,[],options);
uLogL = -uLogL;

% uMLE is the unrestricted maximum likelihood estimate, and uLogL is the loglikelihood maximum.
%
% Impose the restriction to the loglikelihood by setting the corresponding lower and upper bound constraints of \rho to 1. Minimize the negative, restricted loglikelihood.

dof = 1;           % Number of restrictions
rLB = [1 -min(x)]; % Restricted model lower bound
rUB = [1 Inf];     % Restricted model upper bound
[rMLE,rLogL] = fmincon(nLogLGradFun,theta0,[],[],[],[],rLB,rUB,[],options);
rLogL = -rLogL;

% rMLE is the unrestricted maximum likelihood estimate, and rLogL is the loglikelihood maximum.
%
% Use the likelihood ratio test to assess whether the data provide enough evidence to favor the unrestricted model over the restricted model.

[h,pValue,stat] = lratiotest(uLogL,rLogL,dof)

% pValue is close to 0, which indicates that there is strong evidence suggesting that the unrestricted model fits the data better than the restricted model.
% Copyright 2012 The MathWorks, Inc.