% Exercise on study of FC intervals
% fcp 140519
% Compute FC intervals as a function of mean for a N(mu,1) distribution
% Physical values of mu are non-negative
%%
sigma = 1;     % standard deviation of sampling normal
CL = normcdf(1,0,1)-normcdf(-1,0,1);      % "68%" confidence level
alpha = 1-CL;

% for y vs mu plot
mu = 0.01:.001:1;
% for mu vs y plot
%mu = 0.01:.01:3.11;

n = length(mu);
low=-200;
high=0;
ylo = zeros(1,n);
yhi = zeros(1,n);

for i=1:n       % loop over values of mu
    mui=mu(i);
    u_greater = @(log_lambda_alpha)(mui + sqrt(-2*log_lambda_alpha));
    l_greater = @(log_lambda_alpha)(max(0,mui - sqrt(-2*log_lambda_alpha)));
    u_lesser = 0;
    l_lesser = @(log_lambda_alpha)(min(0,mui/2 + log_lambda_alpha/mui));

    PgePos = @(log_lambda_alpha) ...
        (normcdf(u_greater(log_lambda_alpha),mui,sigma) ...
        - normcdf(l_greater(log_lambda_alpha),mui,sigma));
    PgeNeg = @(log_lambda_alpha) ...
        (normcdf(u_lesser,mui,sigma) ...
        - normcdf(l_lesser(log_lambda_alpha),mui,sigma));
    Pdiff = @(log_lambda_alpha) ...
        abs(PgePos(log_lambda_alpha) + PgeNeg(log_lambda_alpha)-CL);
    
    [lcrit, diff, exitflag] = fminbnd(Pdiff,low,high);

    % We have the critical likelihood ratio, now construct the table of
    % observations for which we will accept mu
    yhi(i) = sqrt(-2*lcrit) + mui;
    if l_greater(lcrit) > 0
        ylo(i) = -sqrt(-2*lcrit) + mui;
    else
        ylo(i) = l_lesser(lcrit);
    end
end

%%
% Compute Bayes intervals with uniform prior
xlimits = [-3,3];
y = xlimits(1):.01:xlimits(2);
ny = length(y);
bu = zeros(1,ny);
bl = zeros(1,ny);

for i=1:ny
    Llz = normcdf(0, y(i), sigma);
    CLscaled = CL*(1-Llz); 
    delta = norminv(.5*(CLscaled+1), 0, sigma);
    if (y(i)-delta > 0)
        bu(i) = y(i) + delta;
        bl(i) = y(i) - delta;
    else
        bu(i) = norminv(CLscaled + Llz, y(i), sigma);
        bl(i) = 0.;
    end
end

%%
set(0,'DefaultLineLineWidth',2);
set(0,'DefaultAxesFontSize',15);
set(0,'DefaultTextFontSize',17);

%%
plot( mu,ylo,'b-', mu,yhi,'g-', mu,mu+1,'r--', mu,mu-1,'c--' )
figure;
semilogx( mu,ylo,'b-', mu,yhi,'g-', mu,mu+1,'r--', mu,mu-1,'c--' )
xlim([0,mu(n)]);
%ylim([min(ylo),1.1*max(yhi)+.5]); % min(ylo)=-112.9340  1.1*max(yhi)+.5=2.7
ylim([-10 1.1*max(yhi)+.5])
xlabel('\mu');
ylabel('y_1 and  y_2');

%%
% Now graph the confidence intervals
figure;
yliml = -1;
ylimu = 3;
yc = [yliml,ylimu];
plot(yhi, mu, 'b-', ylo,mu, 'g-', yc-1, yc, 'r--', yc, yc-1, 'c--', y, bl, 'm-', y, bu, 'y-');
xlim(xlimits);
ylim([yliml,ylimu]);
xlabel('y');
ylabel('confidence interval for  \mu');

hold on;
plot(xlimits,[0,0],'color','k');
grid minor
