*! version 6.0.0  15jan1999
program def genhwi, rclass
	version 6
	qui {
        	gettoken AA 0 : 0, parse(" ,")
        	gettoken Aa 0 : 0, parse(" ,")
        	gettoken aa 0 : 0, parse(" ,")

        	confirm integer number `AA'
        	confirm integer number `Aa'
        	confirm integer number `aa'

        	if `AA'<0 | `Aa'<0 | `aa'<0  {
                	di in red "negative numbers invalid"
                	exit 498
        	}

        	syntax [, Label(string) Binvar ]
		tempname obs nA na 
		scalar `obs'= `AA'+`Aa'+`aa'
		scalar `nA'=2*(`AA')+`Aa'
		scalar `na'=2*(`aa')+`Aa'
		local nob=`nA'+1
		if _N < `nA' {
			preserve	
			set obs `nob'
		}
		tempvar nAA nAa naa pobs pcal pval sump
		tempvar fn fnA f2n f2n_na 2x fx
		gen int `nAA'=_n-1
		gen int `nAa'=`nA'-2*(`nAA')
		replace `nAA'=. if `nAa'<0
		#delimit ;
		gen double `pcal'= /*
		*/ lnfact(`obs') - lnfact(2*`obs') + lnfact(`nA') +
		lnfact(2*`obs'-`nA') - lnfact(`nAa')-lnfact((`nA' - `nAa')/2)-
		lnfact(`obs'-(`nA'+`nAa')/2);
		#delimit cr
		replace `pcal'=exp(`pcal'+(`nAa')*ln(2))
		*replace `pcal'=`pcal' * 2^(`nAa')
		gen double `pobs'=`pcal' if `nAa'==`Aa'	
		sort `pobs' `pcal'
		replace `pobs'=`pobs'[_n-1] if `pobs'==. & `pcal'!=.
		gen double  `pval'=`pcal' if `pcal'<=`pobs'
		gen double `sump'=sum(`pval')
		tempname pA pa eAA eAa eaa xAA xAa xaa x2 eobs
		scalar `pA'= (2*(`AA')+`Aa')/(2*`obs')
		scalar `pa'= (2*(`aa')+`Aa')/(2*`obs')
		scalar `eAA'=`obs'*(`pA'^2)	
		scalar `eAa'=`obs'*2*(`pA'*`pa')	
		scalar `eaa'=`obs'*(`pa'^2)	
		scalar `xAA'=((`AA'-`eAA')^2)/`eAA'	
		scalar `xAa'=((`Aa'-`eAa')^2)/`eAa'	
		scalar `xaa'=((`aa'-`eaa')^2)/`eaa'	
		scalar `x2'=`xAA'+`xAa'+`xaa'
		scalar `eobs'=`eAA'+`eAa'+`eaa'
		
		if "`label'"~="" {
			parse "`label'", parse (" ")
			local  naAA= "`1'"
			local  naAa= "`2'"
			local  naaa= "`3'"
		}	
		else {
			local  naAA= "AA"
                	local  naAa= "Aa"
                	local  naaa= "aa"
		}
		local al1=substr("`naAA'",1,1)
		local al2=substr("`naaa'",1,1)
		* noi di in gr _col(6)/*
		*/ "---------- Hardy-Weinberg Equilibrium Test ----------" _n
		/* noi di _n in gr "     Genotype counts" */
		noi di in gr "        Genotype |     Observed        Expected"
		noi di in gr "     ------------+-----------------------------"
		noi di in gr %16s /*
		*/"`naAA'" " |" in ye _col(24) %8.0f `AA' _col(40) %8.2f `eAA' 
		noi di in gr %16s /* 
		*/ "`naAa'" " |" in ye _col(24) %8.0f `Aa' _col(40) %8.2f `eAa' 
		noi di in gr  %16s /*
		*/ "`naaa'" " |" in ye _col(24) %8.0f `aa' _col(40) %8.2f `eaa' 
		noi di in gr "     ------------+-----------------------------"
		noi di in gr _col(12) "total" " |" /*
		*/  in ye _col(24) %8.0f `obs' _col(40) %8.2f `eobs' 
		tempname se1 se2
		if `"`binvar'"'=="" {
			scalar `se1'= /* 
			*/sqrt((1/(2*`obs'))*(`pA'+(`AA'/`obs')-2*(`pA'^2)))
			scalar `se2'= /*
			*/ sqrt((1/(2*`obs'))*(`pa'+(`aa'/`obs')-2*(`pa'^2)))
		}
		else {
			scalar `se1'=sqrt(`pA' * `pa' / (2 * `obs'))
			scalar `se2'=`se1'
			local vartype " (binomial)"
		}
		noi di ""
		/* noi di in gr "     Allele frequencies      " */
		noi di in gr /*
		*/  "          Allele | Observed     Frequency      Std. Err."
		noi di in gr /*
		*/ "     ------------+--------------------------------------"
		noi di in gr /*
		*/ %16s "`al1'" " |" in ye _col(20) %8.0f `nA' /*
		*/ _col(34) %8.4f `pA'  _col(49) %8.4f `se1' "`vartype'"
		noi di in gr /*
		*/ %16s "`al2'" " |" in ye _col(20) %8.0f `na'/*
		*/  _col(34) %8.4f `pa' _col(49) %8.4f `se2' "`vartype'"
		noi di in gr /*
		*/ "     ------------+--------------------------------------"
		noi di in gr _col(12) "total" " |" /*
		*/  in ye _col(20) %8.0f `na'+`nA' _col(34) %8.4f `pA'+`pa' 
		tempname D VD llo lle lr ED
		scalar `D'=(`AA'/`obs')-(`pA'^2)
		scalar `ED'=`D'-(1/(2*`obs'))*(`pA'*(1-`pA')+`D')
		scalar `VD'=(1/`obs')* /*
		*/ ( (`pA'^2) * ((1-`pA')^2) + ((1-2*`pA')^2)*`D' - `D'^2)
		*noi di (`D'-`ED')/sqrt(`VD')
		scalar `llo'= /* 
		*/ (`AA'*ln(`AA'/`obs'))+(`Aa'*ln(`Aa'/`obs'))+ /*
		*/ (`aa'*ln(`aa'/`obs'))
		scalar `lle'= (`AA'*ln(`eAA'/`eobs')) + /*
		*/ (`Aa'*ln(`eAa'/`eobs')) + (`aa'*ln(`eaa'/`eobs'))
		scalar `lr'=2*(`llo'-`lle')
		noi di _n in gr /*
		*/ "     Estimated disequilibrium coefficient (D) = " /*
		*/ in ye %8.4f `D' 
		/*  in gr "  SE= " in ye %5.4f sqrt(`VD') */
		noi di _n in gr "     Hardy-Weinberg Equilibrium Test: "
		noi di in gr /*
		*/ "              Pearson chi2 (" in ye "1" in gr ") = " /* 
		*/ in ye %8.3f `x2' /*
		*/  in gr "  Pr= " in ye %5.4f  chiprob(1,`x2')
		noi di in gr /*
		*/ "     likelihood-ratio chi2 (" in ye "1" in gr ") = " /* 
		*/ in ye %8.3f `lr' /*
		*/  in gr "  Pr= " in ye %5.4f  chiprob(1,`lr')
	
		noi di in gr  /*
		*/ "      Exact significance prob  =               " /*
		*/ in ye %5.4f `sump'[_N]  
		ret scalar  D = `D'
		ret scalar p_exact = `sump'[_N]
		ret scalar lr_p = chiprob(1,`lr')
		ret scalar lr_chi2 = `lr'
		ret scalar chi2_p = chiprob(1,`x2')
		ret scalar df = 1
		ret scalar chi2 = `x2'
}
	end