program define power
    version 4.0

    preserve    
    clear       

    quietly {
        set obs 501
        gen x= .
        gen xbar1=.
        gen xbar2=.
        gen s21=.
        gen s22=.
        gen stat = .
        gen y = .}

    global D_sm12 "   The power of a test is the probability of correctly rejecting H_0, that is, rejecting it when it is false. This lab illustrates this idea. It is exactly the same as the level of significance lab except: 1) All tests in this lab are right-sided, and 2) you must specify how false H_0 is (see the description of D below)."
    wdctl static D_sm12 5 5 300 32

    global D_sm13 "   Choosing D: D is a number between 0 and 100 which indicates how far the true value of the parameter being tested is from the hypothesized value. The value 0 means the true value is in fact the hypothesized value (thus H_0 is true), while the value 100 means the true value is very far from the hypothesized value. The other values of D (25, 50, and 75) fall in between 0 and 100 in how far the true and hypothesized values are apart."
    wdctl static D_sm13 5 40 300 40

    global D_sm14 "Things to Study: The effect of alpha, D, and sample size on the proportion of times H_0 is rejected by the lab. (Power increases as D and/or sample size(s) and/or alpha increase.)"
    wdctl static D_sm14 5 85 300 32


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

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

    global D_sm31 "alpha"
    global D_sm32 ".01 .025 .05 .10"
    wdctl static D_sm31          122 120 20 10
    global D_var4 ".05"
    wdctl ssimple D_var4 D_sm32  122 130 30 50 

    global D_sm4 "D"
    global D_sm5 "0 25 50 75 100"
    wdctl static D_sm4           155 120 20 16 
    global D_var3 "0"
    wdctl ssimple D_var3 D_sm5   155 130 30 60

    global D_sm8 "Statistic"
    global D_sm9 "Z One-sample-t Two-sample-t Chi-square F"
    wdctl static D_sm8            188 120 45 10
    global D_var2 "Z"
    wdctl ssimple D_var2 D_sm9    188 130 55 60

    global D_sm6 "Sample from"
    global D_sm7 "Normal(0,1) Uniform(0,1) Exponential(1)"
    wdctl static D_sm6            245 120 45 10
    global D_var1 "Normal(0,1)"
    wdctl ssimple D_var1 D_sm7    245 130 55 40

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


    global D_b1 "powdr"
    global D_b2 "exit 1234"
    global D_b3 "whelp power"

    cap noi wdlg "Power of a Test" 5 5 310 200

    restore

end

program define powdr
    version 4.0

    local n1=$HJN_V1
    local n2=$HJN_V2
    local alpha=$D_var4
    local D=$D_var3/100.

    if `n1'< 2 | `n1' > 100 {
        sstopbox stop "n1 must be between 2 and 100"
        exit}

    if `n2'< 2 | `n2' > 100 {
        sstopbox stop "n_2 must be between 2 and 100"
        exit}

    local nsamps=1
    if "$D_var2"=="F" | "$D_var2"=="Two-sample-t" {local nsamps=2}

    global D_var91="Normal"
    if "$D_var1"=="Uniform(0,1)" {global D_var91="Uniform"}
    if "$D_var1"=="Exponential(1)" {global D_var91="Exponential"}


    gph open
    gph pen 2

    if "$D_var1"=="Normal(0,1)" { 
        local sig2 = 1.0 
        local mu = 0.0}
    if "$D_var1"=="Uniform(0,1)" { 
        local sig2 = 1./12
        local mu = 0.5}
    if "$D_var1"=="Exponential(1)" { 
        local sig2 = 1.0
        local mu = 1.0}

    local sig = sqrt(`sig2')

    replace xbar1=0
    replace s21=0

    if `nsamps'==1 {
        gph text 1000 16000 0 0 Generating 500 Samples From a $D_var91 Population
        gph text 2000 16000 0 0 and Calculating $D_var2 for Each Sample}
    else {
        gph text 1000 16000 0 0 Generating 500 Pairs of Samples From a $D_var91 Population
        gph text 2000 16000 0 0 and Calculating $D_var2 for Each Pair of Samples}

    gph box 6000 10000 8000 20000 4

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

        local bx2=int(10000 + 10000*(`i'/`n1'))
        gph box 6000 10000 8000 `bx2' 1

        if "$D_var1"=="Uniform(0,1)"     { replace  x = uniform()}
        if "$D_var1"=="Exponential(1)" { replace  x = -log(uniform())}
        if "$D_var1"=="Normal(0,1)"      { replace  x = invnorm(uniform())}

        if "$D_var2"=="Z" | "$D_var2"=="One-sample-t" {
            replace x=x+ `sig'*`D'}

        if "$D_var2"=="Chi-square" {replace x=x*sqrt(1+`D')}

        if "$D_var2"=="F" {replace x=x*sqrt((1+`D')*`sig2')}

        if "$D_var2"=="Two-sample-t" {replace x=x+`D'*`sig'}
        

        replace xbar1=xbar1+x
        replace s21=s21+x^2
        local i=`i'+1
    }

    replace s21=(s21-(xbar1^2/`n1'))/(`n1'-1)
    replace xbar1=xbar1/`n1'

    if `nsamps'==2 {

        replace xbar2=0
        replace s22=0

        gph box 6000 10000 8000 20000 4

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

            local bx2=int(10000 + 10000*(`i'/`n2'))
            gph box 6000 10000 8000 `bx2' 1

            if "$D_var1"=="Uniform(0,1)"     { replace  x = uniform()}
            if "$D_var1"=="Exponential(1)" { replace  x = -log(uniform())}
            if "$D_var1"=="Normal(0,1)"      { replace  x = invnorm(uniform())}

            replace xbar2=xbar2+x
            replace s22=s22+x^2
            local i=`i'+1
        }

        replace s22=(s22-(xbar2^2/`n2'))/(`n2'-1)
        replace xbar2=xbar2/`n2'
    }

    if "$D_var2"=="Z"            {replace stat = (xbar1-`mu')/sqrt(`sig2'/`n1')}
    if "$D_var2"=="One-sample-t" {replace stat = (xbar1-`mu')/sqrt(s21/`n1')}
    if "$D_var2"=="Chi-square"   {replace stat = ((`n1'-1)*s21)/`sig2'}
    if "$D_var2"=="Two-sample-t" {
                    replace s21=((`n1'-1)*s21+(`n2'-1)*s22)/(`n1'+`n2'-2)
                    replace stat = (xbar1-xbar2)/sqrt(s21*((1/`n1')+(1/`n2')))}
    if "$D_var2"=="F"            {replace stat = s21/s22}

    set graphics off
    graph stat, bin(20) xlabel ylabel l1(" ") l2(" ") b1(" ") b2(" ") bbox(3000,0,23000,32000,850,400,0)
    set graphics on

    local ymin = _result(1)
    local ymax = _result(2)        
    local xmin = _result(3)
    local xmax = _result(4)

    replace x = `xmin' + ((_n-1)/500)*(`xmax'-`xmin')
    local ymax=`ymax'/((`xmax'-`xmin')/20)

    replace y = `ymax'*uniform()

    gph close
    gph open
    gph pen 1

    graph y x, s(i) xlab ylab l1(" ") l2(" ") b1(" ") b2(" ") bbox(3000,0,23000,32000,850,400,0)
    gphconv

    global D_var5="Right-Sided"

    local alp = `alpha'
    if "$D_var5"=="Two-Sided" {local alp = `alpha'/2}

    if "$D_var2"=="Z" {
        replace y = exp(-(x^2)/2)/sqrt(2*_pi)
        local critr=invnorm(1.-`alp')
        local critl=-`critr'}

    local ioptt=0
    if "$D_var2"=="One-sample-t" | "$D_var2"=="Two-sample-t" {local ioptt=1}
    if `ioptt'==1 {
        local df=`n1'-1
        if "$D_var2"=="Two-sample-t" {local df=`n1'+`n2'-2}
        pdf_t x y `df'
        local critr=invt(`df',1.-2*`alp')
        local critl=-`critr'}

    if "$D_var2"=="Chi-square" {
        local df=`n1'-1
        pdf_chi2 x y `df'
        local critr=invnchi(`df',0,1.-`alp')
        local critl=invnchi(`df',0,`alp')}

    if "$D_var2"=="F" {
        local df1=`n1'-1
        local df2=`n2'-1
        pdf_f x y `df1' `df2'
        local critr=invfprob(`df1',`df2',`alp')
        local critl=invfprob(`df1',`df2',1.-`alp')}

*    replace y = y*(`xmax'-`xmin')/20
*    summarize y
*    local yy1 = _result(6)
*    if(`yy1' > `ymax') {local ymax=`yy1'}



*    graph stat, bin(20) xlabel ylabel xscale(`xmin',`xmax') bbox(3000,0,23000,32000,850,400,0)/*
*        */ yscale(`ymin',`ymax') l1(" ") l2(" ") b1(" ") b2(" ") shading(4)
*    graph y x, c(l) s(i)  xlabel ylabel xscale(`xmin',`xmax') bbox(3000,0,23000,32000,850,400,0) /*
*        */ yscale(`ymin',`ymax') l1(" ") l2(" ") b1(" ") b2(" ") pen(3)
*    gphconv

    replace xbar1=.
    replace xbar2=.

    if "$D_var5"=="Right-Sided" {
        replace xbar1=x if x>`critr' 
        replace xbar2=y if x>`critr'}

    if "$D_var5"=="Left-Sided" {
        replace xbar1=x if x<`critl' 
        replace xbar2=y if x<`critl'}

    if "$D_var5"=="Two-Sided" {
        replace xbar1=x if x<`critl' | x > `critr' 
        replace xbar2=y if x<`critl' | x > `critr'}

    gph pen 3
    lines x y

    replace s21=0
    segments xbar1 s21 xbar1 xbar2

    gph pen 2
    histgm stat `xmin' `xmax' 20


*    local mcritr=critr
*    local mcritl=critl

    gph pen 3

    local ycr = -`ymax'/6

    if "$D_var5"=="Two-Sided" {
        local mcritr: display %6.3f `critr'
        local mcritl: display %6.3f `critl'
        drtext `critr' `ycr' 1 `mcritr'
        drtext `critl' `ycr' -1 `mcritl'}
    if "$D_var5"=="Right-Sided" {
        local mcrit: display %6.3f `critr'
        drtext `critr' `ycr' 1 `mcrit'}
    if "$D_var5"=="Left-Sided" {
        local mcrit: display %6.3f `critl'
        drtext `critl' `ycr' -1 `mcrit'}


*    drpoint `crit' 0  1000 6
*    drpoint `mcritl' 0 1000 1
*    drpoint `mcritr' 0 1000 1


    if "$D_var5"=="Right-Sided" {
        replace stat=stat>`critr'
        replace stat=sum(stat)}

    if "$D_var5"=="Left-Sided" {
        replace stat=stat<`critl'
        replace stat=sum(stat)}

    if "$D_var5"=="Two-Sided" {
        replace stat=stat>`critr' | stat < `critl'
        replace stat=sum(stat)}


    local ng:  display %5.3f stat[501]/501.
    local ng1: display %6.3f `crit'

    gph pen 2

    if `nsamps'==2 {
        gph text 1000 16000 0 0 Values of $D_var2 for 500 Pairs of Samples (n_1=`n1',n_2=`n2') From a $D_var91 Population, D=$D_var3}
    else {
        gph text 1000 16000 0 0 Values of $D_var2 for 500 Samples (n=`n1') From a $D_var91 Population, D=$D_var3}

    gph text 2000 16000 0 0 Proportion of Times $D_var5 Test With alpha=$D_var4 is Rejected is `ng'

*    if `nsamps'==2 {
*        gph text 3000 16000 0 0 n_1=`n1', n_2=`n2'}
*    else {
*        gph text 3000 16000 0 0 n=`n1'}

    gph close

end

exit