function varargout = computeIntegralDis(disD, probabilityVector)
% compute Monte Carlo integral estimate from distribution
%
% Syntax: value = computeIntegralDis(disD)
%
% Returns the Monte Carlo estimate of the distribution, based on the 
%   function values in disD.valueVector and the probability density in for
%   each sample point in disD.probabilityVector.
%
% [value, absErr, relErr] = computeIntegralDis(..) also returns a measure
%   for the relative and absolute error of the integral value. These are
%   the empirical standard deviation from the single point estimates and
%   the empirical coefficient of variation. Total error bounds are not
%   available and the real error might well be higher.
%
% computeIntegralDis(valueVector, probabilityVecor) is an alternative
%   call to perform the computation.
%
% See also: discretizeDistribution, reduceDistribution

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2013 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU General Public License as published by the 
%   Free Software Foundation, either version 3 of the License, or (at your 
%   option) any later version.
%
% This program is distributed in the hope that it will be useful, but 
%   WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
%   General Public License for more details.
%
% You should have received a copy of the GNU General Public License along 
%   with this program.  If not, see <http://www.gnu.org/licenses/>.

if ~isfield(disD, 'positionMatrix') || ~isfield(disD, 'valueVector') || ...
        ~isfield(disD, 'probabilityVector')
    error('elk:distribution:wrongInput', ['input structure must consist ' ...
        'of the fields positionMatrix, valueVector and probabilityVector']);
end

%% Input handling
if isstruct(disD)
    valueVector = disD.valueVector(:)';
    probabilityVector = disD.probabilityVector(:)';
else
    valueVector = disD(:)';
    probabilityVector = probabilityVector(:)';
end
if numel(valueVector) ~=numel(probabilityVector)
    error('elk:distribution:wrongInput', ['Vector length for values and ' ...
        'probability do not match']);
else
    nSample = numel(valueVector);
end


%% Calculations
% element wise information
estimateVector = valueVector ./ probabilityVector;

% integral value
value = 1/nSample*sum(estimateVector);

% absolute error rating (~68,3% confidence level)
absErr = sqrt(1/nSample/(nSample-1) * sum(...
    (estimateVector - value).^2));

% relative error
relErr = absErr / value;

%% Output handling
varargout{1} = value;
if nargout >= 2
    varargout{2} = absErr;
end
if nargout >= 3
    varargout{3} = relErr;
end
