*! version 1.0.0  14may1995
program define cis
    version 4.0

    preserve                   /*  Save all the current stata stuff */
    clear                      /*  Clear everything so no conflicts */

    global D_sm11 "   This lab shows 1) interpreting confidence intervals, 2) the effect of sample sizes and alpha on the width of intervals, and 3) what happens to our confidence if violations of assumptions are made."
    wdctl static D_sm11 5 5 290 25

    global D_sm12 "   You specify 1) which parameter to study, 2) what sample size(s) to use, 3) how many intervals to draw, 4) what alpha to use, and 5) which population to sample from. Then the lab will repeatedly get random samples and for each one draw the confidence interval."
    wdctl static D_sm12 5 40 290 25

    global D_sm13 "   Suggestions: 1) run the lab several times for the same parameter and population but varying sample size(s) and/or alpha to see how the width of the intervals changes; 2) run the lab for small samples from the Uniform or Exponential and see how the confidence level is different than what it should be."
    wdctl static D_sm13 5 75 290 35

    global D_sm1 "n_1"
    global HJN_V1=20
    wdctl static D_sm1  5 120 12 10 
    wdctl edit   HJN_V1 20 120 16 10 

    global D_sm2 "n_2"
    global HJN_V2=20
    wdctl static D_sm2  40 120 12 10 
    wdctl edit   HJN_V2 55 120 16 10 

    global D_sm31 "alpha"
    global D_sm32 ".01 .025 .05 .10 .25"
    wdctl static D_sm31  105 110 20 10
    global D_var4 ".05"
    wdctl ssimple D_var4 D_sm32 105 120 30 60 

    global D_sm4 "# Ints"
    global D_sm5 "50 100 150"
    wdctl static D_sm4         140 110 20 10 
    global D_var3 "50"
    wdctl ssimple D_var3 D_sm5 140 120 25 48

    global D_sm6 "Parameter"
    global D_sm7 "mu Variance mu_1-mu_2 Var_1/Var_2 pi pi_1-pi_2"
    wdctl static D_sm6           170 110 40 10
    global D_var1 "mu"
    wdctl ssimple D_var1 D_sm7   170 120 55 70

    global D_sm8 "Sample from"
    global D_sm9 "Uniform(0,1) Exponential(1) Normal(0,1) 0-1"
    wdctl static D_sm8           230 110 45 10
    global D_var2 "Normal(0,1)"
    wdctl ssimple D_var2 D_sm9   230 120 55 60

    wdctl button "Run"     5 160 30 14 D_b1
    wdctl button "Close"    40 160 30 14 D_b2
    wdctl button "Help"    75 160 30 14 D_b3 help

    global D_b1 "cisdr"
    global D_b2 "exit 1234"
    global D_b3 "whelp cis"

    cap noi wdlg "Confidence Intervals" $D_dlgx $D_dlgy 300 240

    restore

end

program define cisdr
    version 4.0

/*  Put Stuff from Dialog Boxes into Local Macros: */

    local n1=$HJN_V1
    local n2=$HJN_V2
    local alpha=$D_var4
    local nints=$D_var3

/*  Check Input: */

    if `n1' < 5 | `n1' > 500 | `n2' < 5 | `n2' > 500 {
        sstopbox stop "sample size must be between 5 and 500"
        exit}    

    if `alpha' <= 0. | `alpha' >= 1.0 {
        sstopbox stop "alpha must be between 0 and 1"
        exit}

    if `nints' < 2 | `nints' > 200 {
        sstopbox stop "number of intervals must be between 2 and 200"
        exit}

    if "$D_var1" == "pi" | "$D_var1" == "pi_1-pi_2" {
        if "$D_var2" != "0-1" {
            sstopbox stop "You must sample from a 0-1 population"
            exit}
    }
    else {
        if "$D_var2" == "0-1" {
            sstopbox stop "You cannot sample from a 0-1 population"
            exit}
    }

/*  Get Number of Samples and What to set obs to: */

    local nsamps=2
    local n=max(`n1',`n2')

    if "$D_var1"=="mu" | "$D_var1"=="Variance" | "$D_var1"=="pi" {
        local nsamps=1
        local n=`n1'}

/*  Get critical values to be used in confidence intervals (Z, t, chi-square, F) */

    if "$D_var1"=="mu" { 
        local crit=invt(`n1'-1,1-`alpha')}
    if "$D_var1"=="mu_1-mu_2" {
        local crit=invt(`n1'+`n2'-2,1-`alpha')}
    if "$D_var1"=="pi" | "$D_var1"=="pi_1-pi_2" {
        local crit=invnorm(1-(`alpha'/2))}
    if "$D_var1"=="Variance" {
        local chiu=invnchi(`n1'-1,0,`alpha'/2)
        local chil=invnchi(`n1'-1,0,1-(`alpha'/2))}
    if "$D_var1"=="Var_1/Var_2" {
        local fu=invfprob(`n1'-1,`n2'-1,`alpha'/2)
        local fl=invfprob(`n1'-1,`n2'-1,1-(`alpha'/2))}

/*  Define local macro true = true value: */

    local true=0
    if "$D_var1"=="pi" {local true=0.5}
    if "$D_var1"=="Variance" | "$D_var1"=="Var_1/Var_2" {local true=1}
    if "$D_var1"=="Variance" & "$D_var2"=="Uniform(0,1)" {local true = 1./12.}
    if "$D_var1"=="mu" & "$D_var2"=="Uniform(0,1)" {local true=.5}
    if "$D_var1"=="mu" & "$D_var2"=="Exponential(1)" {local true=1.0}

/*  Specify lower and upper limits of horizontal axis scale: */

    local xlow=0.0  /* These work for "pi" */
    local xup=1.0

    if "$D_var1"=="Variance" & "$D_var2"=="Normal(0,1)" {
        local xlow=0.
        local xup=5.}

    if "$D_var1"=="Variance" & "$D_var2"=="Exponential(1)" {
        local xlow=0.
        local xup=5.}

    if "$D_var1"=="Variance" & "$D_var2"=="Uniform(0,1)" {
        local xlow=0.
        local xup=1.}

    if "$D_var1"=="pi_1-pi_2" {
        local xlow=-1.
        local xup=1.}

    if "$D_var1"=="Var_1/Var_2" & "$D_var2"=="Normal(0,1)" {
        local xlow=0.
        local xup=10.}

    if "$D_var1"=="Var_1/Var_2" & "$D_var2"=="Uniform(0,1)" {
        local xlow=0.
        local xup=10.}

    if "$D_var1"=="Var_1/Var_2" & "$D_var2"=="Exponential(1)" {
        local xlow=0.
        local xup=10.}

    if "$D_var1"=="mu" & "$D_var2"=="Normal(0,1)" {
        local xlow=-4
        local xup=4}
    if "$D_var1"=="mu" & "$D_var2"=="Exponential(1)" {
        local xlow=0
        local xup=5}
    if "$D_var1"=="mu" & "$D_var2"=="Uniform(0,1)" {
        local xlow=0
        local xup=1}

    if "$D_var1"=="mu_1-mu_2" {
        local xlow=-4
        local xup=4}

/*  numin keeps track of intervals containing true value: */

    local numin=0

/*  Open graphics, set obs and generate variables: */

    gph open
    gph pen 1

    quietly {
        set obs 200
        gen x=.
        gen y=.
        gen x1=.
        gen y1=.}


/*  Draw Axes: */

    replace x=`xlow'+(`xup'-`xlow')*uniform()
    replace y=1
    replace y=`nints' in `n'
    graph y x, xlab(`xlow',`true',`xup') ylab(0,`nints') xscale(`xlow',`xup')/*
            */ yscale(0,`nints') s(i) bbox(3000,0,23000,32000,850,400,0)
    gphconv

/*  Draw Vertical Line at True Value: */

    gph pen 3
    replace x1=`true' in 1
    replace x1=`true' in 2
    replace y1=1 in 1
    replace y1=`nints' in 2
    lines x1 y1

/*  Start loop: */

    local percent: display %5.1f 100*(1-`alpha')

*    gph text 500 500 0 -1 `percent'% CI's for $D_var1


    if `nsamps'==1 {
        gph text 500 500 0 -1 `percent'% CI's for $D_var1; Sampling from $D_var2, n=`n1'}

    if `nsamps'==2 {
        gph text 500 500 0 -1 `percent'% CI's for $D_var1; Sampling from $D_var2, n_1=`n1', n_2=`n2'}



    local i=1
    while `i' <= `nints' {

/*    Get first sample and its mean and variance: */

        if "$D_var2"=="Normal(0,1)" {replace x = invnorm(uniform()) in 1/`n1'}
        if "$D_var2"=="Exponential(1)" {replace x = -log(uniform()) in 1/`n1'}
        if "$D_var2"=="Uniform(0,1)" {replace x = uniform() in 1/`n1'}
        if "$D_var2"=="0-1" {replace x=uniform() < .5 in 1/`n1'}

        qui summarize x in 1/`n1'
        local xbar1=_result(3)
        local s21=_result(4)
 
/*    If necessary, get second sample and its mean and variance: */

        if `nsamps'==2 {
            if "$D_var2"=="Normal(0,1)" {replace x = invnorm(uniform()) in 1/`n2'}
            if "$D_var2"=="Exponential(1)" {replace x = -log(uniform()) in 1/`n2'}
            if "$D_var2"=="Uniform(0,1)" {replace x = uniform() in 1/`n2'}
            if "$D_var2"=="0-1" {replace x=uniform() < .5 in 1/`n2'}
        }

        qui summarize x in 1/`n2'
        local xbar2=_result(3)
        local s22=_result(4)

/*  Now get lower and upper limits and put in ll and lu: */

        if "$D_var1"=="mu" {
            local se=sqrt(`s21'/`n1')
            local ll=`xbar1'-`crit'*`se'
            local lu=`xbar1'+`crit'*`se'}

        if "$D_var1"=="pi" {
            local se=sqrt(`xbar1'*(1-`xbar1')/`n1')
            local ll=`xbar1'-`crit'*`se'
            local lu=`xbar1'+`crit'*`se'}

        if "$D_var1"=="Variance" {
            local ll=(`n1'-1)*`s21'/`chil'
            local lu=(`n1'-1)*`s21'/`chiu'}

        if "$D_var1"=="mu_1-mu_2" {
            local s2p=((`n1'-1)*`s21'+(`n2'-1)*`s22')/(`n1'+`n2'-2)
            local se=sqrt(`s2p'*((1/`n1')+(1/`n2')))
            local ll=(`xbar1'-`xbar2')-`crit'*`se'
            local lu=(`xbar1'-`xbar2')+`crit'*`se'}

        if "$D_var1"=="Var_1/Var_2" {
            local sratio=`s21'/`s22'
            local ll=`sratio'/`fu'
            local lu=`sratio'/`fl'}

        if "$D_var1"=="pi_1-pi_2" {
            local se=sqrt((`xbar1'*(1-`xbar1')/`n1')+((`xbar2'*(1-`xbar2')/`n2')))
            local ll=(`xbar1'-`xbar2')-`crit'*`se'
            local lu=(`xbar1'-`xbar2')+`crit'*`se'}

/*  Check for inclusion and set pen number accordingly. Also update inclusion
    counter: */

        gph pen 2
        if `ll' < `true' & `lu' > `true' {
            local numin=`numin'+1
            gph pen 3}

/*  Draw line for this interval. Note we don't allow interval to go outside
    predetermined scale: */

        replace y1=`i' in 1
        replace y1=`i' in 2
        replace x1=max(`ll',`xlow') in 1
        replace x1=min(`lu',`xup') in 2

        lines x1 y1

/*  Keep track of smallest lower limit and largest upper limit (this is a
    temporary thing to help get scales): */

        if `i'==1 {
            local llmin=`ll'
            local lumax=`lu'}
        else {
            if `ll' < `llmin' {local llmin=`ll'}
            if `lu' > `lumax' {local lumax=`lu'}
        }

/*  Increment loop counter and end loop: */

        local i=`i'+1

    }        

/*  Write message at top of screen: */
    
    local llmin: display %6.2f `llmin'
    local lumax: display %6.2f `lumax'

    local numin: display %5.1f 100*(`numin'/`nints')
    gph pen 2

    gph text 2500 1500 0 -1 % of intervals containing $D_var1 = `numin'

/*  Close graphics and drop variables: */

    gph close
    drop y x y1 x1

end
exit