/* THIS PROGRAM COMPUTES A BIAS CORRECTED AND ACCELERATED BOOTSTRAPPED 95% CI FOR A MEAN */.
/* Syntax is BOOTMEAN X where X is a variable in the current data file. Output provides the sample mean, the */.
/* bootstrap estimate of the population mean based on 1000 resamples with replacement, and the lower and upper */.
/* limits of a bias corrected and adjusted 95% confidence interval for the population mean */.
DEFINE bootmean (!POSITIONAL !TOKEN (1)).
set mxloops = 10000.
matrix.
get dd/variables = !1/MISSING = OMIT.
compute reps = 1000.
compute n = nrow(dd).
compute mnb = make(reps,1,0).
compute dat=dd.
compute mni=dd.
/* Here we generate the statistic of interest, MN */.
compute mn=csum(dat)/n.
/* HERE WE GENERATE THE OMITTED CASE ESTIMATES, MNI */.
compute tmp=make(n-1,1,0).
loop #n = 1 to n.
compute b = 0.
loop #m = 1 to n.
do if (#m <> #n).
compute b=b+1.
compute tmp(b,1)=dd(#m,1).
end if.
end loop.
/* Here we generate the statistic of interest */.
compute mni(#n,1)=csum(tmp)/(n-1).
end loop.
/* HERE WE GENERATE THE INFLUENCE STATISTICS, U */.
compute u=(mni-mn)*((n-1)/n).
/* HERE WE GENERATE THE ACCELERATION ESTIMATE, A */.
compute a = csum((u/n)&**3)/(6*((csum((u/n)&**2))&**(3/2))).
/* HERE WE GENERATE THE BOOTSTRAP RESAMPLE ESTIMATES, MNB */.
loop #n = 1 to reps.
loop #m = 1 to n.
compute v=trunc(uniform(1,1)*n)+1.
compute dat(#m,1)=dd(v,1).
end loop.
compute mnb(#n,1)=csum(dat)/n.
end loop.
/* This is the place where the actual statistic of interest is computed */.
compute bmn = csum(mnb)/reps.
/* NOW WE GENERATE THE
compute #pv = mnb < mn.
compute #pv = csum(#pv)/reps.
compute p = #pv.
compute p0=-.322232431088.
compute p1 = -1.
compute p2 = -.342242088547.
compute p3 = -.0204231210245.
compute p4 = -.0000453642210148.
compute q0 = .0993484626060.
compute q1 = .588581570495.
compute q2 = .531103462366.
compute q3 = .103537752850.
compute q4 = .0038560700634.
do if (#pv > .5).
compute p = 1-#pv.
end if.
compute y=sqrt(-2*ln(p)).
compute xp=y+((((y*p4+p3)*y+p2)*y+p1)*y+p0)/((((y*q4+q3)*y+q2)*y+q1)*y+q0).
do if (#pv <= .5).
compute xp = -xp.
end if.
compute zlo = xp-((1.96-xp))/(1+(a*(1.96-xp))).
compute zhi = xp+((xp-(-1.96)))/(1-(a*(xp-(-1.96)))).
compute plo = cdfnorm(zlo).
compute phi = cdfnorm(zhi).
compute lo = trunc(plo*(reps+1)).
compute hi = (reps+1)-trunc((1-phi)*(reps+1)).
do if (lo < 1).
compute lo = 1.
end if.
do if (hi > reps).
compute hi = reps.
end if.
/* NOW WE SORT THE BOOTSTRAPPED ESTIMATE VECTOR */.
compute mnb = {-999;mnb}.
loop #i = 2 to reps+1.
compute ix = mnb(#i,1).
loop #k= #i to 2 by -1.
compute k = #k.
do if (mnb(#k-1,1) > ix).
compute mnb(#k,1)=mnb(#k-1,1).
else if (mnb(#k-1,1) <= ix).
BREAK.
end if.
end loop.
compute mnb(k,1)=ix.
end loop.
compute mnb = mnb(1:reps+1,1).
/* NOW WE FIND THE UPPER AND LOWER LIMITS */.
compute loci = mnb(lo,1).
compute hici = mnb(hi,1).
compute bw={mn, bmn, loci, hici, n}.
print/title = "BIAS CORRECTED AND ACCELERATED BOOTSTRAP MEAN ESTIMATES, 1000 RESAMPLES".
print bw/title = " "/clabels = "Sample" "Bootstrp" "Lo95%CI " "Hi95%CI" "n"/format f9.4.
end matrix.
!END DEFINE.