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