Confidence Intervals for Binomial Proportion (Again): A Quick Note

In Lex’s library of the latest SAS Global Forum 2015 papers, I found an interesting paper by Wu Gong, Jeffreys Interval for One-Sample Proportion with SAS/STAT Software, where SAS MCMC procedure and a so called Random Walk Metropolis Algorithm were implemented to calculate the Jeffreys interval for binomial proportion.

Years ago I wrote several posts on this topic:

Statistical Notes (3): Confidence Intervals for Binomial Proportion Using SAS

Statistical Notes (5): Confidence Intervals for Difference Between Independent Binomial Proportions Using SAS

Confidence Intervals for Binomial Proportion: A SAS 9.4/STAT 12.3 Update

Confidence Intervals for Difference Between Independent Binomial Proportions: A SAS 9.4/STAT 12.3 Update

I’m not a statistician and I might get some time later to dig these new methods in Wu’s paper (I modified his codes little bit to fit my post):

/*0 input data*/
data jeff;
    n=263;
    r=81;
run;

/*1 Random Walk Metropolis*/
%let nmc=1010000;
%let c=0.08;

data PosteriorSample;
    call streaminit(123);
    set jeff;
    p0=0.5;
    retain p0;

    do i = 1 to &nmc.;
        p1=p0+rand(‘normal’,0,&c.);

        do while(p1 lt 0 or p1 gt 1);
            p1=p0+rand(‘normal’,0,&c.);
        end;

        logratio=(1/2-1)*log(p1)+(1/2-1)*log(1-p1)+r*log(p1)+(n-r)*log(1-p1)
        – (1/2-1)*log(p0)-(1/2-1)*log(1-p0)-r*log(p0)-(n-r)*log(1-p0);

        if log(rand(‘uniform’)) <= logratio then do;
            p0=p1;
        end;

        if i gt 10000 and floor(i/20)*20 eq i then do;
            output;
        end;
    end;
    keep i p0;
run;

title "Jeffreys Interval by Random Walk Metropolis Algorithm";
proc univariate data=PosteriorSample noprint;
    var p0;
    output out=PP pctlpts = 2.5 97.5 pctlpre = pct;
run;

proc print data=pp noobs;
    format pct2_5 pct97_5 6.4;
run;
title;

/*2 MCMC*/
title "Jeffreys Credible Interval by PROC MCMC";
ods select PostIntervals;
proc mcmc data=jeff
    seed=123
    outpost=PosteriorSample2
    nbi=10000
    nthin=20
    nmc=1000000
    statistics=(summary interval) diagnostics=none plots=none;
    parms prb 0.5;
    prior prb ~ beta(1/2,1/2);
    model r ~ binomial(n,prb);
run;
title;