*! version 3.0   03/26/92  (STB10: sbe9)
program define briertst
        version 3.0
        local varlist "req min(3) max(3)"
        local if "opt"
        local in "opt"
        local options "Group(integer 10)"
        parse "`*'"
        parse "`varlist'", parse(" ")
        cap assert `1'==0 | `1'==1
        if (_rc) { error 149 }

        tempvar FMD GID DJ INSAMP NGROUP FJ DJ LOG AB DENOM

        quietly {
                sort `2'
                gen byte `INSAMP' = 1 `if' `in'
                replace `INSAMP' = 0 if `1'==. | `2'==. | `3'==.
                gen float `FMD' = sum(`INSAMP')
                local totsamp = `FMD'[_N]
                replace `FMD' = (`FMD'-1)/`FMD'[_N]        /* relative percentile */
                summ `1' if `INSAMP'
                local dbar = _result(3)
                summ `2' if `INSAMP'
                local fbar = _result(3)
                local fvar = _result(4) * (_result(1)-1)/_result(1)

                summ `2' if `INSAMP' & `1'==1
                local f1 = _result(3)
                replace `FMD' = `f1' if `1'==1
                summ `2' if `INSAMP' & `1'==0
                local f0 = _result(3)
                replace `FMD' = `f0' if `1'==0
                summ `FMD' if `INSAMP'
                local sfmin = _result(4) * (_result(1)-1)/_result(1)

                replace `FMD' = (`2'-`1')^2                /* (f(i)-d(i))^2 */
                summ `FMD' if `INSAMP'
                local brier1 = _result(3)

                replace `FMD' = (`3'-`1')^2                /* (f(i)-d(i))^2 */
                summ `FMD' if `INSAMP'
                local brier2 = _result(3)

                local num=`brier1' - `brier2'
                gen AB=(`2' + `3')/2 if `INSAMP'
                qui sum AB if `INSAMP'
                local pi=_result(3)
                local nsq=_result(1)^2
                gen DENOM= ((`2' - `3')^2) * `pi' * (1 - `pi') if `INSAMP'
                replace DENOM=sum(DENOM)
                local VB=(4/`nsq')*DENOM[_N]
                local comp=abs(`num' / sqrt(`VB'))

        }

di in gr "Comparison of two Brier scores" _col(35) in ye %7.4f `comp' /*
 */"; p = "%7.4f 2*(1-normprob(`comp'))
/*

        di `brier'-(`oiv'+`fvar'+`relinla'-2*`sfd1')
        di `brier1'-`sanders'-`relinsm'
        di `sanders'-`oiv'+`murphy'
        di `relinsm'-`fvar'-`relinla'+2*`sfd1'-`murphy' - `gerr'

        assert abs(`brier1'-`sanders'-`relinsm')<.0001
        assert abs(`sanders'-`oiv'+`murphy')<.0001
        assert abs(`relinsm'-`fvar'-`relinla'+2*`sfd1'-`murphy')<.0001

The last assertion does not hold because there is a difference between
brier and brier1 conceptions.  Yates statements hold for the original
Brier score, but the Sanders and Murphy decompositions hold for the
brier1 score.  The difference between brier and brier1 is the
grouping error.  It is probably somewhat more complicated than I
computed here.

*/
end