function [Conf,Conftable,devtable]=Confint(X,group,plots,conflev)
% Construct confidence intervals for each value of a grouping variable
%
%
%  Required input arguments:
%
%            X: numeric vector 
%           group : numeric or logical vector, character array, string
%                   array, or cell array of character vectors, containing
%                   group names. This is the so called grouping variable.
%
%  Optional input arguments:
%
%   plots : Boolean scalar. If true the confidence intervals are shown on
%           the screen. The default value of plots is false.
%  conflev : scalar in the interval (0 1). conflev defines the confidence
%           level which must be used to compute the confidence intervals
%
%  Output:
%
%        Conf:   array of dimension k-by-3 which contains the confidence interval
%               of the means of X for each of the k groups of the grouping
%               variable. First column of Conf refers to the lower values
%               of the confidence intervals while third column refers to
%               the upper values of the confidence intervals. The second
%               column contains the means.
%   Conftable:  table of dimension k-by-3 containing the same information
%               of Conf but with rowNames and colNames
%     devtable: table of dimension 1-by-2. The first element contains the
%               withing groups deviance (in percentage) and the second
%               element the between groups deviance (in percentage). 
%
%
% Examples
%{
    % numeric vector
    strength = [82 86 79 83 84 85 86 87 74 82 ...
                78 75 76 77 79 79 77 78 82 79];
    % grouping variable
    alloy = {'st','st','st','st','st','st','st','st',...
             'al1','al1','al1','al1','al1','al1',...
             'al2','al2','al2','al2','al2','al2'};
        [Conf,Conftable,devtable]=Confint(strength,alloy)
%}


%{
    % In this case the confidence intervals are shown on the screen
    % and 99 per cent confidence intervals are used

    % numeric vector
    strength = [82 86 79 83 84 85 86 87 74 82 ...
                78 75 76 77 79 79 77 78 82 79];
    % grouping variable
    alloy = {'st','st','st','st','st','st','st','st',...
             'al1','al1','al1','al1','al1','al1',...
             'al2','al2','al2','al2','al2','al2'};
    % Confint is called with four elements
        [Conf,Conftable,devtable]=Confint(strength,alloy,true,0.99)
%}

%% Beginning of code

% Controllo sugli argomenti di input
if ~isnumeric(X)
    error('MiaFunzione:Confint:WronInt','Attenzione il vettore di input X deve essere numerico')
end

if nargin==2
    plots=false;
    conflev=0.95;
elseif nargin==3
    conflev=0.95;
end

% Controllo che conflev sia uno scalare compreso tra 0 e 1
if ~isscalar(conflev)
    error('MiaFunzione:Confint:Wronconflev','Attenzione conflev deve essere uno scalare')
end
if conflev>=1 || conflev <=0
    error('MiaFunzione:Confint:Wronconflev','Attenzione conflev deve essere un numero compreso tra 0 e 1')
end

% Nel caso in cui l'utente passi un valore [] per plots
if isempty(plots)
    plots=false;
end

% unigroup=cell che contiene le etichette dei diversi gruppi
unigroup=unique(group);

% inizializzazione di devNEI e devFRA
devNEI=0;
devFRA=0;
ngroups=length(unigroup);
% conflev = scalare associato all'ampiezza dell'intervallo di confidenza
zalpha=norminv(1-(1-conflev)/2);
Mglobal=mean(X);

% inizializzazione della matrice ConfInt
% Prima colonna di Confint = estremo inferiore intervallo di confidenza del
% gruppo j
% Seconda colonna di Confint = Media del gruppo j
% Terza colonna di Confint = estremo superiore intervallo di confidenza del
% gruppo j
Conf=zeros(ngroups,3);

for j=1:ngroups
    booj=strcmp(group,unigroup(j));
    % nj = numerosità del gruppo j
    nj=sum(booj);
    % Xj= vettore di lunghezza nj che contiene le unità appartenenti al
    % gruppo j
    Xj=X(booj);
    % varj = varianza (corretta) delle unità appartenenti al gruppo j
    varj=var(Xj);
    % stdj= standard deviation (s corretto) delle unità appartenenti ai
    % gruppo j
    stdj=sqrt(varj);
    
    devNEI=devNEI+varj*(nj-1);
    % Mj=media del gruppo j
    Mj=mean(Xj);
    devFRA=devFRA+(Mj-Mglobal)^2 *nj;
    
    % Costruzione dell'intervallo di confidenza per le unità appartenenti
    % al gruppo j
    Conf(j,:)=[Mj-zalpha*stdj/sqrt(nj) Mj  Mj+zalpha*stdj/sqrt(nj)];
end

% devTOT = devianza totale
devTOT=devNEI+devFRA;


% costruzione degli output devtable e Conftable
devtable=array2table(100*[devNEI devFRA]/devTOT,'VariableNames',{'DevNei' 'DevTra'});
Conftable=array2table(Conf,'RowNames',unigroup,'VariableNames',{'Low' 'Mean' 'Up'});

% Se plots è true allora gli intervalli di confidenza vengono rappresentati
% graficamente utilizzando la funzione errorbar
if plots==true
    errorbar((1:ngroups)',Conf(:,2),(Conf(:,3)-Conf(:,2)),'o')
    % Le etichette contenente i valori della variabile classificatoria
    % vengono aggiunti nella parte alta del grafico
    text([0.2; 0.5; 0.8],1.05*ones(ngroups,1),unigroup,'Units','normalized','HorizontalAlignment','Center')
    xlim([0.5 ngroups+0.5])
    set(gca,'XTickLabel','')
    ylabel(['Confidence interval with ' num2str(100*conflev) '% level'])
end

end