function [Occ, VND, sigma2, VS, fitted_values, residuals] = occreml(VT)
% occreml - Restricted maximum likelihood of neuroreceptor occupancies
%
% Usage:
%
%   occreml(VT)
%
% Arguments:
%
%   VT - matrix of total volumes of distribution (VT). Each row must have
%        the VT values of a region of interest (ROI). The 1st column must
%        have the baseline VT values, the 2nd column the 1st-postdose VT
%        values, the 3rd column the 2nd-postdose VT values, etcetera.
%
% Value:
%
%   Occ: Vector of postdose neuroreceptor occupancies.
%   VND: Non-displaceable volume of distribution.
%   sigma2: Variance.
%   VS: Vector of ROIs specific volumes of distribution.
%   fitted_values: Matrix of fitted VT.
%   residuals: Matrix of VT residuals.
%
% References:
%
%   Radua J, Bullich S, Lopez N and Catafau AM. Restricted maximum
%   likelihood estimation of PET neuroreceptor occupancy in the absence
%   of a reference region. Medical Physics, 2011, 38:2558-2562.
%
% Example:
%
%   % Input data: occreml_example %
%   occreml_example = [0.72 0.35 0.47 ; 0.96 0.43 0.62 ; 0.75 0.34 0.50 ; 0.44 0.26 0.31 ; 0.39 0.28 0.30]
%   [Occ, VND, sigma2, VS, fitted_values, residuals] = occreml(occreml_example)

    % Get parameters %
    [n, k] = size(VT);
    k = k - 1;
    if k < 1
        stop('No postdose VT found');
    end

    % Minimize the minus log likelihood %
    function [result] = mll(PAR)

        % Prepare parameters %
        Occ = PAR(1:k);
        VND = PAR(k + 1);
        sigma2 = PAR(k + 2) .^ 2;
        VS = PAR((k + 3):(k + 3 + n - 1));

        % Check restrictions %
        if min(Occ) < 0 || max(Occ) > 1 || min(VS) < 0
            result = 9e+30;
            return;
        end

        % Prepare other %
        CO = [1, 1 - Occ]';
        mVT = mean(VT);
        mVT2 = mean(VT .^ 2);

        % Return the mll %
        result = (k + 1) / 2 * n * log(2 * pi * sigma2) + ...
            1 / sigma2 * ( ...
                (k + 1) / 2 * n * VND .^ 2 + ...
                1 / 2 * sum(CO .^ 2) * sum(VS .^ 2) + ...
                sum(CO) * VND * sum(VS) - ...
                VS * VT * CO - ...
                n * VND * sum(mVT) + ...
                n / 2 * sum(mVT2) ...
            );
    end
    reml = fminunc(@mll, [ones(1,k) * 0.5, 0, std(VT(:,1)), VT(:,1)'], optimset('LargeScale','off'));

    % Prepare return %
    Occ = reml(1:k);
    CO = [1, 1 - Occ];
    VND = reml(k + 1);
    sigma2 = reml(k + 2) .^ 2;
    VS = reml((k + 3):(k + 3 + n - 1))';
    fitted_values = VND + VS * CO;
    residuals = VT - fitted_values;

end