DEFINE !BINCINT (N !TOKENS(1) / X !TOKENS(1) /CONF !TOKENS(1) ).
*.
* .
* Given data from a binomial experiment with X = number of
* successes, and N = number of trials, this program returns
* the exact (not asymptotic) upper and lower confidence limits .
* .
title 'Exact Confidence Limits for a Binomial Parameter'.
* .
* Sample Usage: !BINCINT N=10 X=3 CONF=95.0 .
*.
* Mike Lacy, Dec 1996.
* Create a dummy case for SPSS to work with .
new file .
input program.
compute dummy = 1 .
end case.
end file.
end input program.
* .
***********************************************.
* Constants, user modifiable .
compute toler = 1.0e-4. /* desired relative precision .
compute maxiter = 40. /* Maximum iterations allowed below .
* .
string lconverg (a3) uconverg (a3) . /* flags to indicate convergence
* .
* Pick up macro parameters and compute with them as necessary .
compute N = !N .
compute x = !X .
compute alpha2 = (1 - (!CONF/100))/2 .
* Some minor computations.
compute pyes = x/n.
compute se = sqrt(pyes * (1-pyes)/N).
*.
*.
*******************************************.
* Find limits, with first pass doing lower limit
and second pass doing upper limit .
* .
loop dolow = 1 to 0 by -1.
do if dolow .
* Initialize for lower limit .
+ compute lowbound = 0.0 . /*anything too low will work
+ compute upbound = pyes . /*anything too high will work
* Generate two guesses for the lower limit.
+ compute OldGuess = pyes.
+ compute guess = pyes - probit(alpha2) * se.
+ if (guess >= upbound) or (guess <= lowbound)
guess = (lowbound + upbound)/2.0.
* Determine probabilities for initial guesses.
+ compute p = 1 - cdf.binom(x-1,n,guess).
+ compute OldP = 1- cdf.binom(x-1,n,OldGuess).
+ compute lconverg = 'no'.
*.
else. /* Upper Limit .
*.
* Initialize for upper limit .
+ compute lowbound = pyes . /*anything too low will work
+ compute upbound = 1.0 . /*anything too high will work
* Generate two guesses for the upper limit.
+ compute OldGuess = pyes.
+ compute guess = pyes + probit(1-alpha2) * se.
+ if (guess >= upbound) or (guess <= lowbound)
guess = (lowbound + upbound)/2.0.
* Determine probabilities for initial guesses.
+ compute p = cdf.binom(x,n,guess).
+ compute OldP = cdf.binom(x,n,OldGuess).
+ compute uconverg = 'no'.
* .
end if. /* End initialization .
*.
*.
*. Go ahead and find the limit, checking first for extreme cases.
*.
do if (x = 0) and dolow . /* Observed X at lower limit, nothing to do
+ compute lcl = 0.0.
+ compute lconverg = 'yes'.
else if (x=N) and not dolow. /* Observed X at upper limit, nothing to do
+ compute ucl = 1.0 .
+ compute uconverg = 'yes' .
else .
* Use iteration to find the limit The technique is the secant .
* method (empirical derivative), supplemented with bisection .
*.
+ compute numiter = 0 .
+ compute diff = p - alpha2 .
*.
* ITERATION LOOP .
+ loop .
+ compute numiter = numiter + 1 .
* Reset bounds on iteration: different for upper and lower limit .
+ do if (diff < 0) and dolow .
+ compute lowbound = guess.
+ else if (diff <0) and not dolow.
+ compute upbound = guess .
+ else if (diff >= 0) and dolow .
+ compute upbound = guess .
+ else if (diff >= 0) and not dolow .
+ compute lowbound = guess .
+ end if .
* .
* Construct next guess and probability, while protecting
* against bad guesses and wild secant slopes .
+ compute slope = (OldP - P)/(OldGuess - Guess) .
+ compute OldP = P .
+ compute OldGuess = Guess .
+ do if ((Abs(slope) < 1e-5 ) or (Abs(slope) > 1e+5)) .
+ compute guess = (upbound + lowbound)/2.0 .
+ else .
+ compute guess = guess - diff/slope .
+ if (guess >= upbound) or (guess <= lowbound)
guess = (lowbound + upbound)/2.0.
+ end if .
+ do if dolow.
+ compute p = 1- cdf.Binom(x-1,n,Guess) .
+ else if not dolow .
+ compute p = cdf.Binom(x,n,Guess) .
+ end if.
+ compute diff = p - alpha2 .
*.
+ end loop if ( (abs(diff) <= toler * alpha2)
or (numiter > maxiter) ) .
* .
end if . /* if had to find the limit by iteration .
*.
do if dolow .
+ compute lcl = guess .
+ compute l_iter = numiter .
+ if (numiter < maxiter) lconverg = 'yes'.
else if not dolow .
+ compute ucl = guess .
+ compute u_iter = numiter .
+ if (numiter < maxiter) uconverg = 'yes'.
end if . /* if dolow .
end loop . /* loop dolow = 1 to 0 for lower and upper limit .
*.
format pyes (f7.5) ucl (f7.5) lcl (f7.5)
x (f6.0) n (f6.0) l_iter (f6.0) u_iter (f6.0).
list x n lcl pyes ucl l_iter lconverg u_iter uconverg.
!ENDDEFINE .
* .
!BINCINT N=10 X=3 CONF=95.0 .
* .
* X = observed # of "Yes" .
* N = number of trials .
* lcl = lower confidence limit .
* pyes = observed proportion of "Yes" .
* ucl = upper confidence limit .
* l_iter = number of iterations used to find lower limit .
* l_converg = Did iteration for lower limit converge? .
* 'u_iter = number of iterations used to find upper limit .
* 'u_converg = Did iteration for upper limit converge? .
* .
* .
exe .