## Copyright (C) 2012 Martin Šíra
##
## 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 2 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/>.

## -*- texinfo -*-
## @deftypefn {Function File} {[@var{sci}, @var{lq}, @var{rq}] =} scovint (@var{data}, @var{cp}, [ @var{verbose} ] )
## Calculates shortest coverage interval for given coverage 
## probability cp from vector of samples of distribution. 
## The interval can be asymmetric. Thus the output is uncertainty
## according to Evaluation of measurement data - Supplement 1 
## to the "Guide to the expression of uncertainty in measurement"
## - Propagation of distributions using a Monte Carlo method, JCGM, 
## 2008. In the case of multiple shortest coverage intervals, the one 
## symmetric to expectation (mean) is selected.
##
## The function finds shortest coverage interval only such covering the
## expectation. This means for multipeaked PDF there can exist shorter
## coverage interval, however GUM requires only unimodal (single peaked)
## PDF.
##
## Input variables:
## @table @samp
## @item @var{data} - vector of numbers
## @item @var{cp} - coverage interval probability
## @item @var{verbose} - displays output and figures
## @end table
##
## Output variables:
## @table @samp
## @item @var{sci} - shortest coverage interval
## @item @var{lq} - value of output quantity at the left side of shortest coverage interval
## @item @var{rq} - value of output quantity at the right side of shortest coverage interval
## @end table
##
## Example:
## generate random numbers according to normal distribution and calculate 
## uncertainty of the result (compare it with standard deviation times two):
## @example
## a=randn(1,1e5);
## [sci, lq, rq] = scovint(a, 0.6827, 1)
## 2.*std(a)
## @end example
## @end deftypefn

## Authors: Martin Šíra <msiraATcmi.cz>, Stanislav Mašláň
## Created: July 2012
## Version: 1.12
## Keywords: coverage interval GUM

function [sci, lq, rq] = scovint(data, cp, verbose=0)

        % --------------------------- check input parameters %<<<1

        % check No of arguments
        if ( nargin < 2 )
                print_usage();
        endif
        % check if input data is vector 
        if ( isvector(data)~=1 )
                error ("scovint: DATA must be vector");
        endif
        % check cp
        if ( !isscalar(cp) )
                if ( !isnumeric(cp) )
                        error ("scovint: CP must be a scalar number");
                endif
        endif
        % check correct cp
        if ( cp<=0 || cp>=1 )
                error ("scovint: CP must be greater than 0 and smaller than 1");
        endif
        % check correct verbose
        if ( !isscalar(verbose) )
                error ("scovint: verbose must be scalar");
        endif

        % --------------------------- inner constants %<<<1

        % output quantity resolution (defines resolution of result):
        % XXX POSSIBLY CHANGE THIS VALUES IN THE CASE OF IMPROPER RESULTS:
        % XXX BUT THE CALCULATION WILL TAKE LONGER TIME
        cdfxlength = 500;
        probvalraster = 0.0001;

        % --------------------------- preliminary values %<<<1

        % data length
        N = length(data);

        % reshape data to vertical
        data = reshape(data,N,1);

        % --------------------------- cumulative distribution function: %<<<1

        % see figure 1 of 1st GUM supplement

        % output quantity vector (x axis of the cdf plot):
        cdfx(:,1) = linspace(min(data),max(data),cdfxlength);
        % probability vector (y axis of the cdf plot) normalized:
        cdfy = interp1(sort(data)',[0:N-1]',cdfx(1:end))/(N-1);

        % --------------------------- coverage interval length: %<<<1

        % see figure 7 of 1st GUM supplement
        % expectation:
        expe = mean(data);
        % probability of expectation:
        p_expe = interp1(cdfx,cdfy,expe);

        % left hand probability vector (in interval <0,p_expe>):
        % (NOTE! some multi peaked PDF can generate shortest coverage interval
        % not covering expectation, however GUM supplement 2 requires only
        % single peaked PDF, therefore it is safe to use p_expe as a limit
        % for left hand probability vector)
        lhp = [0:probvalraster:p_expe];
        % generate right hand probability vector for every lhp:
        rhp = cp + lhp;
        % limit rhp to interval <p_expe,0>
        id = rhp > p_expe & rhp < 1;          
        rhp = rhp(id);
        % limit also lhp according rhp:
        lhp = lhp(id);

        % output quantity for every lhp and rhp:
        lhq = interp1(cdfy, cdfx, lhp, 'extrap');
        rhq = interp1(cdfy, cdfx, rhp, 'extrap');

        % calculate vector of lengths of coverage intervals:
        loci = rhq - lhq;                   

        % --------------------------- select (almost) shortest lengths of coverage intervals: %<<<1

        % Shortest length of coverage interval ideally should be only one. However
        % lengths of coverage intervals forms a noisy curve therefore a limit is set 
        % and all lengths of coverage intervals below this limit are considered as
        % shortest lengths of coverage intervals.
        % This problem is best seen on rectangular distribution - all lengths of shortest
        % intervals should be equal, however the really shortest one can be very asymmetric 
        % to the expectation and this is not desired.

        % get limit below which all lengths of coverage intervals are considered as shortest:
        % (span of loci multiplied by 1e-5 and add 10 rounding errors for the safety)
        limit = (max(loci) - min(loci))*1e-5 + 10*eps*mean(loci);
        % get indexes of selected shortest lengths of coverage intervals
        selid = find(loci <= min(loci) + limit);

        % get shortest lengths of coverage intervals:
        selloci = loci(selid);                               
        % quantities matching shortest lengths of coverage intervals:
        sellhq = lhq(selid);
        selrhq = rhq(selid);
        % probabilities matching shortest lengths of coverage intervals:
        sellhp = lhp(selid);
        selrhp = rhp(selid);

        % possible shortest interval lengths
        selsloci = selrhq - sellhq;

        % --------------------------- find most symetric shortest length of coverage interval: %<<<1

        % in the case of multiple shortest lengths of coverage intervals, find most symetric results:
        tmp = abs(abs(sellhq - expe) - abs(selrhq - expe));
        symsel = min(tmp);
        symselid = find(tmp == symsel);
        % if more symetric results, use middle one:
        if(length(symselid) > 1)
          symselid = symselid(round(length(symsel)/2));
        endif

        % get all final values:
        symsellhq = sellhq(symselid);
        symselrhq = selrhq(symselid);
        symsellhp = sellhp(symselid);
        symselrhp = selrhp(symselid);
        symselloci = selsloci(symselid);

        % output result:
        sci = symselloci;
        lq = symsellhq;
        rq = symselrhq;


        % --------------------------- plotting: %<<<1

        % plot of the cumulative distribution function:
        if verbose
                figure
                % cdf:
                plot(cdfx,cdfy,'-x')
                hold on
                xl=xlim;
                yl=ylim;
                % show expectation:
                plot([expe expe], [yl(1) yl(2)], '-g');
                % lines to show selected shortest coverage interval:
                % left hand lines - from probability to quantity:
                plot([xl(1) symsellhq], [symsellhp symsellhp], '-r')
                plot([symsellhq symsellhq], [symsellhp yl(1) ], '-r');
                % right hand lines - from probability to quantity:
                plot([xl(1) symselrhq], [symselrhp symselrhp], '-r');
                plot([symselrhq symselrhq], [symselrhp yl(1)], '-r');
                legend('cdf','expectation','shortest cov. int.');
                xlabel('output quantity')
                ylabel('probability')
                title('cumulative distribution function');
                hold off
        endif

        % plot of the length of coverage interval
        if verbose
                figure
                hold on
                plot(lhp,loci,'o-')
                plot(sellhp, selloci, 'xg', 'linewidth',4,'markersize',10)
                plot(symsellhp, symselloci, 'xr', 'linewidth',4, 'markersize',10)
                legend('length of coverage interval', 'shortest l. of c. i.', 'selected s. l. of c. i.'); 
                title(['length of ' num2str(cp) ' coverage interval'])
                xlabel('left hand probability')
                ylabel('length of coverage interval')
                hold off
        endif

endfunction

% --------------------------- tests: %<<<1

%!shared a, b, c, sci, lq, rq
%! a=[-10:0.1:10];
%! b=stdnormal_pdf(a);
%! c=[];
%! for i=1:length(b)
%!         c=[c ones(1,  round(b(i).*1000)  ).*a(i)];
%! endfor
%! [sci, lq, rq] = scovint(c, 0.6827,0);
%!assert(sci, std(c).*2, 0.005);
%!assert(lq<0);
%!assert(rq>0);
%!assert(sci, rq-lq, 2*eps);
%%% test is random here!:
%! c=gamrnd(1,2,1,1e5);
%! [sci, lq, rq] = scovint(c, 0.6827,0);
%!assert(sci, gaminv(0.6827,1,2), 0.02);
%!assert(lq,0,0.01);
%!assert(rq, sci+lq, 0.02);

% --------------------------- demo: %<<<1
%!demo
%! % chosen coverage interval probability:
%! cp=0.6827;
%! % generate points at which distribution will be evaluated:
%! a=[0:0.01:10];
%! % generate probability density function of gamma distribution:
%! b=gampdf(a,2,2);
%! % prepare variable:
%! c=[];
%! % generate samples of the gamma distribution:
%! for i=1:length(b)
%! % for every point of pdf b(i) generate appropriate number of values:
%!         c=[c ones(1,  round(b(i).*1000)  ).*a(i)];
%! endfor
%! % calculate shortest coverage interval:
%! [sci, lq, rq] = scovint(c, cp, 1);
%! % calculate twice the standard deviation:
%! stddouble=std(c).*2;
%! % plot histogram figure:
%! figure
%! hist(c,100);
%! hold on
%! % get figure y axis limits:
%! yl=ylim;
%! % calculate expectation (mean):
%! expe=mean(c);
%! % calculate standard deviation:
%! stdev=std(c);
%! % plot the rest:
%! plot([expe expe], [yl(1) yl(2)], 'g--', 'linewidth',3);
%! plot([expe-stdev expe-stdev], [yl(1) yl(2)], '-y', 'linewidth',3);
%! plot([lq lq], [yl(1) yl(2)], '-r', 'linewidth',3);
%! plot([rq rq], [yl(1) yl(2)], '-r', 'linewidth',3);
%! plot([expe+stdev expe+stdev], [yl(1) yl(2)], '-y', 'linewidth',3);
%! legend('histogram','mean','standard deviation','shortest cov. interv.');
%! xlabel('output quantity');
%! ylabel('count');
%! title('gamma probability density function, A=2, B=2');
%! hold off
%! % 
%! % RESULT: see on probability density function figure - yellow standard
%! % deviation interval is symmetric around expectation (mean), but red
%! % shortest coverage interval is assymetrical and shorter!
%! % RESULT: compare shortest length of coverage interval (sci)
%! % with twice the standard deviation (stddouble):
%! sci
%! stddouble

% vim modeline: vim: foldmarker=%<<<,%>>> fdm=marker fen ft=octave textwidth=1000
