program define nbin
    version 4.0

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

    quietly{
        set obs 101
        gen x1=.
        gen x11=.
        gen x12=.
        gen x2=.
        gen x21=.
        gen x22=.
        gen y1=.
        gen y2=.
        gen y3=.
        gen zero=.}

    global D_sm11 "   If we sample (with replacement) from a 0-1 population (having proportion of 1s denoted by pi), then the number of times we must draw elements from the populations until we get the nth one is denoted by X and is called a Negative Binomial Random Variable and its distribution is called the Negative Binomial Distribution."
    wdctl static D_sm11 5 5 290 32

    global D_sm12 "   In this lab, you choose n and pi and the labs draws a bargraph of the distribution."
    wdctl static D_sm12 5 40 290 8

    global D_sm13 "An Example: Let X be the number of times we flip a coin until we get the 5th head. This would be Negative Binomial with n=5 and pi=.5."
    wdctl static D_sm13 5 60 290 16

    global D_sm14 "Notes: 1) The average number of draws needed to get the nth one is (n-1) + n*(1-pi)/pi. 2) You can vary n and pi to see how these parameters affect the distribution."
    wdctl static D_sm14 75 120 200 32

    global D_sm6 "n"
    global D_sm7 "1 2 3 4 5 10"
    wdctl static D_sm6         5 100 10 10
    global D_var1 "5"
    wdctl ssimple D_var1 D_sm7 5 110 25 70

    global D_sm8 "pi"
    global D_sm9 ".20 .30 .40 .50 .60 .70 .80"
    wdctl static D_sm8         40 100 20 10
    global D_var2 ".50"
    wdctl ssimple D_var2 D_sm9 40 110 25 70

    wdctl button "Run"     115 160 30 14 D_b1
    wdctl button "Close"   150 160 30 14 D_b2
    wdctl button "Help"    185 160 30 14 D_b3 

    global D_b1 "nbindr"
    global D_b2 "exit 1234"
    global D_b3 "whelp nbin"

    nbindr

    cap noi wdlg "The Negative Binomial Distribution" $D_dlgx $D_dlgy 300 240

    restore
end

program define nbindr
    version 4.0

    gph open
    gph pen 1

    local n      = $D_var1
    local pi     = $D_var2
    local mu     = `n'+ `n' * (1.-`pi')/`pi'
    local var    = `n'*(1-`pi')/((`pi')^2)

    local xmin   = max(`n', int(`mu'-3*sqrt(`var')))
    local xmax   = min(int(`mu' + 4 * sqrt(`var'))+1,70)


*    gph text 1000 16000 0 0 `n' `pi' `mu' `var' `xmin' `xmax'

    replace x1=.
    replace x2=.
    replace y1=.
    replace y2=.
    replace y3=.

*    gph text 2000 16000 0 0 `N' `D' `n' `p' `nobs'

    replace x1 = _n in `xmin'/`xmax'
    replace y1 = exp(lngamma(x1)-lngamma(`n')-lngamma(x1-`n'+1)+`n'*log(`pi')+(x1-`n')*log(1-`pi'))

    summarize y1
    local ymax = _result(6)
*    gph text 3000 16000 0 0 `y1max'

    local nm1=`xmin'-1
    local np1=`xmax'+1
    graph y1 x1, s(i) xlab(`nm1',`np1') ylab  l1(" ") l2(" ") b1(" ") b2(" ") /*
            */ bbox(3500,0,23000,32000,850,400,0)
    gphconv

    local yy1=-`ymax'/10.
    local yy2=1.5*`yy1'
    gph pen 2
    local i=`xmin'
    while `i' <= `xmax' {
        if 2*int(`i'/2)==`i' {drtext `i' `yy1' 0 `i'}
        else {drtext `i' `yy2' 0 `i'}
        local i = `i'+1}

    replace x11=x1-.5
    replace x12=x1+.5
    replace zero=0

    gph pen 1
    segments x11 zero x11 y1
    segments x12 zero x12 y1
    segments x11 y1 x12 y1
    segments x11 zero x12 zero

    gph text 1000 16000 0 0 The Negative Binomial Distribution (n=`n', pi=`pi')

    gph close

end
exit