program define dsamp
    version 4.0

    preserve
    clear
    quietly{
        set obs 500
        gen x=.
        gen p=.
    }

    local r1=60
    local c1=5
    local c2=25

    global D_sm1 " This lab plots a probability histogram, simulates data from the distribution, and then superimposes a histogram of the simulated data on the probability histogram. You specify k = the number of possible values for the distribution and numbers w1, ..., wk which are weights specifying the likelihood of the possible values (the probability of a value is its weight divided by the sum of all the weights)."
    wdctl static D_sm1 5 5 240 50

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

        if `i'==7 {local r1=60}

        if `i' > 6 {
            local c1=45
            local c2=65}

        global D_w`i'  "w`i'"
        global w`i'=1

        wdctl static D_w`i'   `c1' `r1' 15 8
        wdctl edit   w`i'     `c2' `r1' 18 8

        local r1=`r1'+10
        local i=`i'+1

    }

    global C_1 "k"
    global C_3 "2 3 4 5 6 7 8 9 10 11 12"
    wdctl static C_1          140 60 10 10
    global C_2 "2"
    wdctl ssimple C_2 C_3     150 60 20 100

    global N_1 "n"
    global N_3 "50 100 150 200 250 300 350 400 450 500"
    wdctl static N_1          100 60 10 10
    global N_2 "100"
    wdctl ssimple N_2 N_3     110 60 20 100

    global DB_rad "Special Cases"
*    wdctl static DB_rad 5 120 60 50 blackframe
*    wdctl static DB_rad 20 120 50 10 center

    wdctl static DB_rad 5 120 60 10
    global DB_chk=1
    wdctl radbegin "Equally likely" 10 130 60 10 DB_chk
    wdctl radio    "Fair-dice"      10 140 60 10 DB_chk
    wdctl radio    "Fair-die"       10 150 60 10 DB_chk
    wdctl radend   "User-specified" 10 160 60 10 DB_chk

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

    global D_b1 "dsampdr"
    global D_b2 "exit 1234"
    global D_b3 "whelp dsamp"

    cap noi wdlg "Sampling From Discrete Distributions" 5 5 260 200

    restore

end

program define dsampdr
    version 4.0

    local k=$C_2
    local n=$N_2



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

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

    if $DB_chk==1 {
        local i=1
        while `i' <= 12 {
            global w`i'=1
            local i=`i'+1
        }
    }

    if $DB_chk==2 {
        local k=12
        global C_2=12
        global w1=0
        global w2=1
        global w3=2
        global w4=3
        global w5=4
        global w6=5
        global w7=6
        global w8=5
        global w9=4
        global w10=3
        global w11=2
        global w12=1}

    if $DB_chk==3 {
        local k=6
        global C_2=6
        global w1=1
        global w2=1
        global w3=1
        global w4=1
        global w5=1
        global w6=1}

    local tot=0

    local i=1
    while `i' <= `k' {
        quietly replace p = ${w`i'} in `i'
        if p[`i'] < 0 {
            sstopbox stop "weights cannot be negative"
            exit}
        local tot=`tot'+p[`i']
        local i=`i'+1
    }

    quietly replace p=p/`tot'


*    discsam `k' `n' p x

    tempvar x1 x2 x3 zero expect
    local kp1=`k'+1
    gen `x1'=_n in 1/`kp1'
    gen `x2'=.
    gen `x3'=.
    gen `expect'=. in 1/`k'
    gen `zero'=0 in 1/`k'

    gph open
    summarize p in 1/`k'
    local pmax=1.25*_result(6)
    local pmin=_result(5)
    replace p=0 in `kp1'

    graph p `x1' in 1/`kp1', s(i) xlab(0,`kp1') ylab yscale(0,`pmax') l1(" ") l2(" ") b1(" ") b2(" ") bbox(5000,0,23000,32000,850,400,0)
    gphconv

    replace p=. in `kp1'
    replace `x1'=. in `kp1'
    local yy1=-`pmax'/10

    gph pen 2

    local i=1
    while `i' <= `k' {
        drtext `i' `yy1' 0 `i'
        local i=`i'+1}

*   Draw wide bars:

    replace `x1'=`x1'-.5 in 1/`k'
    replace `x2'=`x1'+1 in 1/`k'

    gph pen 1
    segments `x1' `zero' `x1' p
    segments `x2' `zero' `x2' p
    segments `x1' p `x2' p
    segments `x1' `zero' `x2' `zero'

    local pm=`pmax'-`yy1'/2
    local i=1
    while `i' <= `k' {
        local pi: display %3.2f p[`i']
        local xi=`x1'[`i']+.5
        local xi=int($AX*`xi'+$BX)
        gph text 3000 `xi' 0 0 `pi'
        local pi=`n'*p[`i']
        replace `expect'=`pi' in `i'
        local pi: display %3.1f `pi'
        gph text 3700 `xi' 0 0 `pi'
*        drtext `xi' `pm' 0 `pi'
        local i=`i'+1
    }

    gph pen 3
    replace x=.
    replace x=uniform() in 1/`n'
    replace p=sum(p)
    count if x<p[1] 
    replace `x3'= _result(1) in 1
    local i = 2
    while `i' <= `k' {

        local im1=`i'-1
        count if x<=p[`i'] & x > p[`im1']
        replace `x3' = _result(1) in `i'
        local i=`i'+1}

    replace `expect'=(`x3'-`expect')^2/`expect'
    replace `expect'=sum(`expect')
    local ee=`expect'[`k']
    local df=`k'-1
    local pv: display %5.2f chiprob(`df',`ee')
    local ee: display %5.2f `ee'


    gph pen 2
    gph text 1000 16000 0 0 Theoretical (wide bars) and Observed (narrow bars) Histograms 
    gph text 2000 16000 0 0 for a Discrete Probability Distribution (Chi-sq = `ee', pval = `pv')
    gph pen 3

    replace p=`x3'/`n' in 1/`k'

    replace `x1'=`x1'+.25
    replace `x2'=`x2'-.25
    segments `x1' `zero' `x1' p
    segments `x2' `zero' `x2' p
    segments `x1' p `x2' p
    segments `x1' `zero' `x2' `zero'

    local pm=`pm'+`yy1'/2
    local i=1
    while `i' <= `k' {
        local pi: display %3.2f p[`i']
        local xi=`x1'[`i']+.25
        local xi=int($AX*`xi'+$BX)
        gph text 5100 `xi' 0 0 `pi'
        local pi=`x3'[`i']
        gph text 4400 `xi' 0 0 `pi'
*        drtext `xi' `pm' 0 `pi'
        local i=`i'+1
    }
    wpupdate
end