bca.m

这是一个三面阵列的数学工具软件

function [lo,up]=BCa(Bsamp,BB,BBB,CIlevel);
% function [lo,up]=BCa(Bsamp,BB,BBB,CIlevel);;
% Input:
%       Bsamp: contains vector with sample solution
%       BB: contains bootstrap parameters (one column per bootstrap)
%       BBB: contains expanded data parameters (one column per bootstrap)
%       CIlevel: level of confidence interval, e.g., .95, or .90
% Output:
%       lo and up are cvectors with BCa lower and upper bounds
%
% uses norminv.m and normcdf.m from statistics toolbox

[m,n_boot]=size(BB);
mm=size(Bsamp,1);
if m~=mm,disp(' pbootstraps do not have same size as sample');return;end;
[mm,n]=size(BBB);
if m~=mm,disp(' parameter matrices do not have same numbers of rows');return;end;

% compute bias parameters
W=BB<Bsamp*ones(1,n_boot);
z=norminv(sum(W')'/n_boot);

% compute acceleration parameters
num=sum((BBB-Bsamp*ones(1,n))'.^3);
denom=6*(sum((BBB-Bsamp*ones(1,n))'.^2)).^1.5;
a=num'./denom';

% determine BCa upper and lower bounds
zCIlo=norminv((1-CIlevel)/2);
zCIup=-zCIlo;
a_lo=round(n_boot*normcdf(z+(z+zCIlo)./(1-a.*(z+zCIlo))));
a_up=round(n_boot*normcdf(z+(z+zCIup)./(1-a.*(z+zCIup))));
w=a_lo==0;  % correction for cases where first percentile happens to become number 0
a_lo(w)=1;
for j=1:m
    BBsort=sort(BB')';
    lo(j)=BBsort(j,a_lo(j));
    up(j)=BBsort(j,a_up(j));
end;