% Exercise on study of FC intervals
% fcp 140502
%
% Look at coverage of deciding on 1 or 2 sided intervals based on data.
% Here we do it "analytically".
%%
set(0,'DefaultLineLineWidth',2);
set(0,'DefaultAxesFontSize',14);
set(0,'DefaultTextFontSize',18);

%%
sigma = 1;     % standard deviation of sampling normal
CL68 = normcdf(1,0,1)-normcdf(-1,0,1);      % "68%" confidence level
CL1 = CL68;     % confidence level for upper limit
%CL1=0.9;
CL2 = CL68;     % confidence level for two-sided interval
nsigma1 = norminv(CL1,0,1);     % number of sigma added to sampling
                                % for upper limit
u = nsigma1*sigma;
nsigma2 = norminv(1-(1-CL2)/2,0,1); % number of sigma added and subtracted 
                                    % for two-sided interval
t = nsigma2*sigma;                                   
pvalue = 0.1;   % p-value for deciding what kind of interval to quote
%pvalue = 0.001;
ycrit = norminv(1-pvalue,0,sigma);  % critical value for upper limit
                                    % vs two-sided interval
mu = 0:.01:5;
pul=normcdf(ycrit,mu,sigma);

%%
plot(mu,pul,'b', mu,1-pul,'g');
xlabel('Parameter value')
ylabel('Probability of interval type')
grid on

%%
n = length(mu);
covul = zeros(1,n);
covts = zeros(1,n);
cdfycrit = pul;
cdfmuminusu = normcdf(mu-u,mu,sigma);
cdfmupt = normcdf(mu+t,mu,sigma);
cdfmumt = normcdf(mu-t,mu,sigma);

%%
i = ycrit > mu -u;
covul(i) = cdfycrit(i)-cdfmuminusu(i);
i = ycrit < mu + t & ycrit > mu - t;
covts(i) = cdfmupt(i)-cdfycrit(i);
j = ycrit < mu - t;
covts(j) = cdfmupt(j)-cdfmumt(j);

%%
figure;
plot(mu,covul,'b', mu,covts,'g', mu,covul+covts,'r');
xlabel('Parameter value');
ylabel('Covering probability');
ylim([0 1]);
grid on
