function [f,xi,bw,fn] = histDENSITY(dat,xi,bw,ke)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%HISTDENSITY  generate a scaled probability density function
% Function generates a pdf, but scales it so the values match that obtained by a histogram with the same binsize
% This is useful if you want to plot it alongside a histogram or if you want the units to be 'stuff' instead
% of probability.
%
% USAGE:
%         f = histDENSITY(dat,xi,bw,ke)
%
% INPUT:
%         dat - data points, vector or matrix, will be converted to column vector
%         xi - vector of points at which to estimate the pdf
%         bw - (optional) binwidth to use
%         ke - (optional) kernel to use, see 'help fitdist' to see a list of kernels, default is normal kernel
%
% OUTPUT:
%    f - scaled pdf
%    xi - the locations where the pdf was estimated (same as input if an input was given)
%    bw - binwidth used (same as input if an input was given) 
%    fn - sum normalised histogram
%
% EXAMPLES:
%
% See also: FITDIST KSDENSITY

% HISTORY:
% version 1.0.0, Release 19/02/19 Initial release
%
% Author: Roddy Grieves
% UCL, 26 Bedford Way
% eMail: r.grieves@ucl.ac.uk
% Copyright 2018 Roddy Grieves

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INPUT ARGUMENTS CHECK
% deal with input variables
inps = {'dat','xi','bw','ke'};
vals = {'0','linspace(min(dat(:)),max(dat(:)),10)','0','''normal'''};
reqd = [1 0 0 0];
for ff = 1:length(inps)
    if ~exist(inps{ff},'var')
        if reqd(ff)
            error('ERROR: vargin %s missing... exiting',inps{ff});            
        end        
        eval([inps{ff} '=' vals{ff} ';']);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FUNCTION BODY
dat = dat(:);
xi = xi(:);

% get PDF
if ~bw
    pd = fitdist(dat,'kernel','kernel',ke);
    bw = pd.BandWidth;
else
    pd = fitdist(dat,'kernel','kernel',ke,'width',bw);   
    bw = pd.BandWidth;    
end

% generate and scale to match histogram
n = sum(~isnan(dat));
binwidth = xi(2)-xi(1); % find the width of each bin
area = n * binwidth;
f = area * pdf(pd,xi);

fn = f(:) ./ numel(dat(:));