function STAT = stats(THIS,THETA,varargin)
% stats  Evaluate selected statistics of an MCMC chain.
%
% Syntax
% =======
%
%     S = stats(POS,THETA,...)
%     S = stats(POS,THETA,LOGPOST,...)
%
% Input arguments
% ================
%
% * `POS` [ poster ] - Posterior simulator, `poster`, object used to
% generate the `THETA` chain.
%
% * `THETA` [ numeric ] - MCMC chain generated by the `arwm` function.
%
% * `LOGPOST` [ numeric ] - Vector of log posterior densities generated by
% the `arwm` function; `LOGPOST` is not necessary if you do not request
% `'mdd'`, the marginal data density.
%
% Output arguments
% =================
%
% * `S` [ struct ] - Struct with the statistics requested by the user.
%
% Options
% ========
%
% * `'estTime='` [ `true` | *`false`* ] - Display and update the estimated time
% to go in the command window.
%
% * `'histBins='` [ numeric | *10* ] - Number of histogram bins if `'hist'` is
% in the `'output'` option.
%
% * `'hpdiCover='` [ numeric | *90* ] - Requested percent coverage of the
% highest probability density interval if `'hpdi'` is in the `'output'`
% option.u
%
% * `'mddGrid='` [ numeric | *0.1:0.1:0.9* ] - Points between 0 and 1 over
% which the marginal data density estimates will be averaged, see
% Geweke (1999).
%
% * `'output='` [ char | cellstr ] - List of output statistics.
%
% * `'progress='` [ `true` | *`false`* ] - Display progress bar in the command
% window.
%
% * `'prctile='` [ numeric | *[10,90]* ] - Requested percentiles if
% `'prctile'` is in the `'output'` option.
%
% Description
% ============
%
% List of available statistics:
%
% --Overall statistics
%
% * `'cov='` -- Covariance matrix of the parameter estimates.
% * `'mdd='` -- Minus the log marginal data density.
%
% --Parameter-specific statistics
%
% * `'chain='` - The entire simulated chain of parameter values
% * `'mean='` - Average.
% * `'std='` - Std deviation.
% * `'median='` - Median.
% * `'mode='` - Mode based on the histogram function outputs.
% * `'hpdi='` - Highest probability density interval with coverage given by
% the option `'hpdi'`.
% * `'hist='` - Histogram bins and counts, with the number of bins `'histbins'`.
% * `'ksdensity='` - The x- and y-axis points for kernel-smoothed posterior
% density.
% * `'bounds='` - Lower and upper bounds set up by the user.
% * `'prctile='` - Percentiles given by the option `'prctile'`.
%
% Example
% ========
%

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

if ~isempty(varargin) && isnumeric(varargin{1})
    logpost = varargin{1};
    varargin(1) = [];
else
    logpost = [];
end

options = passvalopt('poster.statistics',varargin{:});

if ischar(options.output)
    options.output = strtrim(regexp(options.output,'\w+','match'));
end

list = {'chain','cov','mean','median','mode','mdd','std', ...
    'hpdi','hist','bounds','ksdensity','prctile'};

if ischar(options.output)
    options.output = regexp(options.output,'\w+','match');
end

output = struct();
for i = 1 : length(list)
    output.(list{i}) = any(strcmpi(list{i},options.output));
end

isfile = ischar(THETA);

if output.mdd && isempty(logpost) && ~isfile
    utils.error('poster', ...
        ['Vector of log posterior densities must be entered ', ...
        'if lmdd is requested.']);
end

%**************************************************************************

STAT = struct();

if isfile
    [p,t] = fileparts(THETA);
    inputfile = fullfile(p,t);
    master = struct([]);
    mat = {};
    npar = NaN;
    N = NaN;
    dochkposterfiles();
else
    [npar,N] = size(THETA);
end

if output.mean || output.cov || output.std || output.mdd
    thetamean = nan(npar,1);
    domean();
end

if options.progress
    progress = progressbar('IRIS poster.arwm progress');
elseif options.esttime
    eta = esttime('IRIS poster.arwm is running');
end

for i = 1 : npar
    name = THIS.paramList{i};
    
    if isfile
        theta = nan(1,master.DRAW);
        dogettheta(i);
    else
        theta = THETA(i,:);
    end
    
    if output.mode || output.hist
        [histcount,histbins] = hist(theta,options.histbins);
    end
    
    if output.chain
        STAT.chain.(name) = theta;
    end
    if output.mean
        STAT.mean.(name) = thetamean(i);
    end
    if output.median
        STAT.median.(name) = median(theta);
    end
    if output.mode
        pos = histcount == max(histcount);
        STAT.mode.(name) = histbins(pos);
    end
    if output.std
        STAT.std.(name) = ...
            sqrt(sum((theta - thetamean(i)).^2) / (N-1));
    end
    if output.hpdi
        [low,high] = tseries.myhpdi(theta,options.hpdicover,2);
        STAT.hpdi.(name) = [low,high];
    end
    if output.hist
        STAT.hist.(name) = {histbins,histcount};
    end
    if output.prctile
        STAT.prctile.(name) = prctile(theta,[options.prctile],2);
    end
    if output.bounds
        STAT.bounds.(name) = [THIS.lowerBounds(i),THIS.upperBounds(i)];
    end
    if output.ksdensity
        [x,y] = doksdensity();
        STAT.ksdensity.(name) = [x,y];
    end
    
    if options.progress
        update(progress,i/npar);
    elseif options.esttime
        update(eta,i/npar);
    end
    
end

% We subtract the mean from `THETA`. The original `THETA` is not available
% any longer after this point.
if output.cov || output.mdd
    sigma = nan(npar);
    docov();
end

if output.cov
    STAT.cov = sigma;
end

if output.mdd
    uuu = nan(1,N);
    douuu();
    STAT.mdd = domdd();
end

% Nested functions.

%**************************************************************************
    function domean()
        if isfile
            for ii = 1 : npar
                theta = nan(1,N);
                dogettheta(ii);
                thetamean(ii) = sum(theta) / N;
            end
        else
            thetamean = sum(THETA,2) / N;
        end
    end
% domean().

%**************************************************************************
    function docov()
        if isfile
            sigma = zeros(npar);
            for ii = 1 : master.SAVECOUNT
                thistheta = mat{ii}.THETA;
                for jj = 1 : npar
                    thistheta(jj,:) = thistheta(jj,:) - thetamean(jj);
                end
                sigma = sigma + thistheta * thistheta.' / N;
            end
        else
            for ii = 1 : npar
                THETA(ii,:) = THETA(ii,:) - thetamean(ii);
            end
            sigma = THETA * THETA.' / N;
        end
    end
% docov().

%**************************************************************************
    function d = domdd()
        % mdd  Modified harmonic mean estimator of minus the log marginal data
        % density; Geweke (1999).
        
        % Copyright (c) 2010-2012 Troy Matheson & Jaromir Benes.
        logdetsigma = log(det(sigma));
        
        % Compute g(theta) := f(theta) / post(theta) for all thetas,
        % where f(theta) is given by (4.3.2) in Geweke (1999).
        if isfile
            logpost = nan(1,N);
            dogetlogpost();
        end
        logg = -(npar*log(2*pi) + logdetsigma + uuu)/2 - logpost;
        
        % Normalise the values of the g function by its average so that the
        % later sums does not grow too high. We're adding `avglogg` back
        % again.
        avglogg = sum(logg) / N;
        logg = logg - avglogg;
        
        d = [];
        for pr = options.mddgrid(:).'
            crit = chi2inv(pr,npar);
            index = crit >= uuu;
            if any(index)
                tmp = sum(exp(-log(pr) + logg(index))) / N;
                d(end+1) = log(tmp) + avglogg; %#ok<AGROW>
            end
        end
        d = -mean(d);
    end
% domdd().

%**************************************************************************
    function douuu()
        invsigma = inv(sigma);
        if isfile
            pos = 0;
            for ii = 1 : master.SAVECOUNT
                thistheta = mat{ii}.THETA;
                for jj = 1 : npar
                    thistheta(jj,:) = thistheta(jj,:) - thetamean(jj);
                end
                for jj = 1 : size(thistheta,2)
                    pos = pos + 1;
                    uuu(pos) = ...
                        thistheta(:,jj).' * invsigma * thistheta(:,jj); %#ok<MINV>
                end
            end
        else
            % `THETA` is already demeand at this point.
            for jj = 1 : N
                uuu(jj) = THETA(:,jj).' * invsigma * THETA(:,jj); %#ok<MINV>
            end
        end
    end
% douuu().

%**************************************************************************
    function [x,y] = doksdensity()
        lb = THIS.lowerBounds(i);
        ub = THIS.upperBounds(i);
        if isinf(ub) && ~isinf(lb) && lb ~= 0
            if lb > 0
                lb = 0;
            else
                ub = 1e10;
            end
        elseif isinf(lb) && ~isinf(ub) && ub ~= 0
            if ub < 0
                ub = 0;
            else
                lb = -1e10;
            end
        end
        try
            [y,x] = ksdensity(theta,'support',[lb,ub]);
        catch %#ok<CTCH>
            [~,y,x] = thirdparty.kde(theta,2^10,lb,ub);
        end
    end
% doksdensity().

%**************************************************************************
    function dochkposterfiles()
        master = load(inputfile);
        valid = isstruct(master) ...
            && isfield(master,'SAVECOUNT') ...
            && isnumericscalar(master.SAVECOUNT) ...
            && isfield(master,'DRAW') ...
            && isnumericscalar(master.DRAW) ...
            && isfield(master,'PLIST') ...
            && iscellstr(master.PLIST);
        valid = valid && isequal(THIS.paramList,master.PLIST);
        N = master.DRAW;
        npar = length(THIS.paramList);
        if valid
            ndraw = 0;
            for ii = 1 : master.SAVECOUNT
                thisfile = sprintf('%s%g',inputfile,ii);
                mat{ii} = matfile(thisfile);
                thiswho = who('-file',thisfile);
                valid = valid ...
                    && length(thiswho) == 2 ...
                    && any(strcmp(thiswho,'THETA')) ...
                    && any(strcmp(thiswho,'LOGPOST'));
                [npar,ntheta] = size(mat{ii},'THETA');
                [~,nlogpost] = size(mat{ii},'LOGPOST');
                valid = valid ...
                    && ntheta == nlogpost ...
                    && npar == length(master.PLIST);
                ndraw = ndraw + ntheta;
            end
            valid = valid && ndraw == master.DRAW;
        end
        if ~valid
            utils.error('poster', ...
                'The posterior simulation file ''%s'' is invalid.');
        end
    end
% dochkposterfiles().

%**************************************************************************
    function dogettheta(i)
        count = 0;
        for ii = 1 : master.SAVECOUNT
            [~,n] = size(mat{ii},'THETA');
            theta(1,count+(1:n)) = mat{ii}.THETA(i,:);
            count = count + n;
        end
    end
% dogettheta().

%**************************************************************************
    function dogetlogpost()
        count = 0;
        for ii = 1 : master.SAVECOUNT
            [~,n] = size(mat{ii},'LOGPOST');
            logpost(1,count+(1:n)) = mat{ii}.LOGPOST;
            count = count + n;
        end
    end
% dogetlogpost().

end
