st: RE: Simplification of mata syntax to determine CI possible?

[Nick Cox <n.j.cox@durham.ac.uk>]

I'd recommend reading the paper below before doing any work in this territory. The alternatives mentioned, which have long been implemented in -ci-, work well. I can't see any need to loop to get an individual binomial c.i. In any case, stepping uniformly thorough

   pp = range(0.000001,0.999999,1/99999)

sounds a dodgy thing to do. 

Interval Estimation for a Binomial Proportion
Lawrence D. Brown, T. Tony Cai, and Anirban DasGupta
Source: Statist. Sci. Volume 16, Issue 2 (2001), 101-133. 
Abstract
We revisit the problem of interval estimation of a binomial proportion. The erratic behavior of the coverage probability of the standard Wald confidence interval has previously been remarked on in the literature (Blyth and Still, Agresti and Coull, Santner and others). We begin by showing that the chaotic coverage properties of the Wald interval are far more persistent than is appreciated. Furthermore, common textbook prescriptions regarding its safety are misleading and defective in several respects and cannot be trusted. 

This leads us to consideration of alternative intervals. A number of natural alternatives are presented, each with its motivation and context. Each interval is examined for its coverage probability and its length. Based on this analysis, we recommend the Wilson interval or the equal-tailed Jeffreys prior interval for small n and the interval suggested in Agresti and Coull for larger n. We also provide an additional frequentist justification for use of the Jeffreys interval

Nick 
n.j.cox@durham.ac.uk 


-----Original Message-----
From: owner-statalist@hsphsun2.harvard.edu [mailto:owner-statalist@hsphsun2.harvard.edu] On Behalf Of Dirk Enzmann
Sent: 06 July 2011 11:58
To: statalist@hsphsun2.harvard.edu
Subject: st: Simplification of mata syntax to determine CI possible?

I would like to know whether it is possible to simplify the syntax / 
computations of the following mata syntax:

The syntax calculates likelihood based confidence intervals of a 
binomial parameter given sample size n, successes x and confidence level 
ci (I hope the lines will not get broken in the mail):

* ---- start ---------
mata:
void propci(real scalar x, real scalar n, real scalar ci)
{
   real rowvector pp
   real scalar k
   real rowvector lik
   real scalar ci_l
   real scalar ci_u
   real scalar i
   pp = range(0.000001,0.999999,1/99999)
   k = exp(1)^(-invchi2(1, ci/100)/2)
   lik = binomialp(n,x,pp)
   lik = lik/max(lik)
   ci_l = .
   ci_u = 0
   for (i=1; i <= 100000; i++) {
     if (lik[i,1] > k & ci_l == .) {
       ci_l = pp[i,1]
     }
     else if (lik[i,1] > k & ci_u < pp[i,1]) {
       ci_u = pp[i,1]
     }
     else if (ci_u > 0 & lik[i,1] < k) {
       i = 100000
     }
   }
   st_numscalar("r(ci_l)", ci_l)
   st_numscalar("r(ci_u)", ci_u)
}
end
* ---- end -----------

The 13 lines between -lik = lik/max(lik)- and -st_numscalar("r(ci_l)", 
ci_l)- (exclusive) can be (re)written in R using only one line of commands:

# --- start R snippet -------

conf.int=range(pp[lik > k])

# --- end R snippet ------

I wonder whether something similar is possible using mata.

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/