matlab/bca.m

   1 function [lower,upper]=bca(estimate,bootstat,jackstat,low,up)\r
   2 % Beräknar BCA-konfidensintervall enligt Efron\r
   3 % Anrop [lower,upper]=bca(estimate,bootstat,jackstat,low,up)\r
   4 % Parametrar:\r
   5 % estimate=originalestimat\r
   6 % bootstat=matris av bootstrapestimat\r
   7 % jackstat=matris av jackknifeestimat\r
   8 % low=undre konfidensgrad i procent\r
   9 % up=övre konfidensgrad i procent\r
  10 % Output:\r
  11 % lower=undre konfidensgräns\r
  12 % upper=övre konfidensgräns\r
  13 \r
  14 % Beräkning av hjälpstorheten zhat \r
  15 [a,b]=size(bootstat);\r
  16 antal=sum(bootstat<=(ones(a,1)*estimate))\r
  17 zhat=norminv(antal/a);\r
  18 % Beräkning av hjälpstorheten ahat\r
  19 %[nrow,ncol]=size(x);\r
  20 %if nrow<ncol\r
  21 %       x=x';\r
  22 %       nrow=ncol ;\r
  23 %end\r
  24 %% Jackknife estimates in s (Lasse)\r
  25 %for j=1:nrow\r
  26 %       if j==1\r
  27 %               values=[2:nrow];\r
  28 %       elseif j==nrow\r
  29 %               values=[1:nrow-1];\r
  30 %       else    \r
  31 %               values=[1:j-1, j+1:nrow];\r
  32 %       end\r
  33 %       s(j,:)=feval(fun,x(values,:));\r
  34 %end\r
  35 s=jackstat;\r
  36 mn=mean(s);                             % theta_hat_dot (Lasse)\r
  37 [b c]=size(s);\r
  38 d=s-ones(b,1)*mn;\r
  39 sum2=sum(d.^2);                         % - täljare\r
  40 sum3=sum(d.^3);                         % del av nämnare\r
  41 ahat=-sum3./(6*(sum2.^(1.5)));\r
  42 % Beräkning av korrigerade konfidensgrader enligt Efron \r
  43 n=norminv(low/100);\r
  44 alfa1=100*normcdf(zhat+(zhat+n)./(1-ahat.*(zhat+n)));\r
  45 n=norminv(up/100);\r
  46 alfa2=100*normcdf(zhat+(zhat+n)./(1-ahat.*(zhat+n)));\r
  47 % Beräkning av konfidensgränser med percentilmetoden\r
  48 lower=diag(prctile(bootstat,alfa1))';\r
  49 upper=diag(prctile(bootstat,alfa2))';\r