Chapter 5: computing confidence intervals

To understand the basic character of the shortest 95% confidence range (interval), we will develop a ‘brute-force’ method. This method will essentially step through each point in the posterior distribution, and from there compute the length along the abscissa needed to encompass 95% of the probability mass. The shortest of these is our approximation of the shortest 95% confidence interval. The approximation is limited by increment size.

mu=10; sd=12;

mus=linspace(-60,60,1001)';

post=npdf(mus,mu,sd);

Ninterp=201;

%%% step through the posterior

cs=cumtrapz(mus,post);

post=post./cs(end);

cs=cumtrapz(mus,post); %normalize the posterior

itest=0; inow=0; iend=zeros(size(cs)); c0=zeros(size(cs)); %initialize index

while inow<length(mus), cnow=c0;

cnow(inow+1:end)=cumtrapz(mus(inow+1:end),post(inow+1:end));

itest=find(cnow>=.95,1,'first')-1;

if and(itest<length(mus),~isempty(itest)),

inow=inow+1;

iend(inow)=itest;

else break; end, end

startend=[[1:inow]' iend(1:inow)];

ranges=diff(startend,1,2);

imins=find(ranges==min(ranges)); Ushaped=0;

for nmin=1:length(imins),

iminmax=startend(imins(nmin),1); iminmax(1,2)=iminmax(1)+1;

vmins=linspace(mus(iminmax(1)),mus(iminmax(2)),Ninterp+2);

vmins=vmins(2:end-1)';

vmaxs=nan(Ninterp,1);

pmins=interp1(mus,post,vmins);

for inow=1:Ninterp,

ctrap=cumtrapz([vmins(inow); mus(imins(2):end)],[pmins(inow); post(imins(2):end)]);

vmaxs(inow,1)=interp1(ctrap,[vmins(inow); mus(imins(2):end)],0.95); end

lens=vmaxs-vmins;

indmin=find(lens==min(lens));

if ~or(indmin==1,indmin==length(lens)),

Ushaped=1;

break; end, end

if Ushaped, ['success'], else ['error: no min found'], end

CIapprox=[vmins(indmin) vmaxs(indmin)]; CI=[-13.52 33.52]; %icdf('normal',[.025 .975],mu,sd);

figure(1); clf;

subplot(2,1,1); plot(vmins,lens,'ko:','MarkerSize',8);

subplot(2,1,2); hold on; plot(mus,post/max(post),'k.');

for n=1:2, plot(CIapprox(n)*[1 1],[0 1],'b-','LineWidth',2); plot(CI(n)*[1 1],[0 1],'k:','LineWidth',2); end