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 .