I'm attempting to reproduce in Stata a technique from a statistical coding textbook ("Analysis of Categorical Data with R", Bilder and Loughlin).

I ultimately want to show visually how near confidence intervals for a binomial distribution (using the Wald technique) approach the idealized 95% for a range of probability values - 0.001(0.0005)0.999). In other words, something like the graphic attached below (which was generated in R using the author's code).

I'm really trying to iterate the sample code below across the range of 0.001(0.0005)0.999, whereas the below code only gets me pi=0.001 and I would like to get it for the whole range.

clear
set obs 1998
gen w = rbinomial(40,0.001)
gen pi_hat = w/40
gen var_wald = (pi_hat*(1-pi_hat)/40)
gen lower = pi_hat + (invnormal(.025) * sqrt(var_wald))
gen upper = pi_hat - (invnormal(.025) * sqrt(var_wald))
gen pmf=0
replace pmf=binomialp(40,w,pi)
gen save = 1 if pi>lower & pi<upper

gen CI = save/1998

The iterative techniques I've tried encompass many variations of the following, but I invariably get endless loops or broken code that doesn't give me what I want.

forvalues j = 1/1998{
forvalues i = 0.001(0.0005)0.999{
replace pi = `i'
replace w = rbinomial(40,pi)
replace pi_hat = w/40
replace var_wald = (pi_hat*(1-pi_hat)/40)
replace lower = pi_hat + (invnormal(.025) * sqrt(var_wald))
replace upper = pi_hat - (invnormal(.025) * sqrt(var_wald))
replace pmf=0
replace pmf=binomialp(40,w,pi)
replace save = 0
replace save = 1 if pi>lower & pi<upper
}
replace CI = save/1998 in `j'
}

My goal is for the variable CI to contain a list of all the unique proportions describing how often the Wald confidence interval contained the "true" pi, for a given value of pi.

This is my first time posting in the group, although I've followed for years to learn the program better. I do apologize if I've made mistakes in formulating this question.

Thanks in advance,


Jack McDonnell
Array