version 12


capture program drop vx_dbel

program vx_dbel, rclass
	version 12
	syntax varlist  [, alpha(real 0.05) repsnum(real 10000) sided(string) method(string)]

	preserve

	quietly putmata datmatx=(`varlist'), replace

	mata: ncols=cols(datmatx)
	drop _all
	getmata (ncols)=ncols, replace force
	local ncolsflag=ncols

	//set default for string
	if (strmatch(strlower("`method'"),"")) 					local method="interpolate"
	if (strmatch(strlower("`sided'"),"") & `ncolsflag'==1 ) local sided="one.sample"
	if (strmatch(strlower("`sided'"),"") & `ncolsflag'==2 ) local sided="two.sample"
	if (strmatch(strlower("`sided'"),"") & `ncolsflag'==3 ) local sided="three.sample"

	//facilitate change to requested syntax of only two.sided
	if (strmatch(strlower("`sided'"),"two.sided") & `ncolsflag'==1 ) local sided="one.sample"
	if (strmatch(strlower("`sided'"),"two.sided") & `ncolsflag'==2 ) local sided="two.sample"
	if (strmatch(strlower("`sided'"),"two.sided") & `ncolsflag'==3 ) local sided="three.sample"

	//facilitate change to requested syntax by addition of less
	if (strmatch(strlower("`sided'"),"less") & `ncolsflag'==1 ) mata: datmatx=datmatx:*-1
	if (strmatch(strlower("`sided'"),"less") & `ncolsflag'==2 ) mata: datmatx=datmatx:*-1

	//display "`varlist'"
	//display "`method'"
	//display "`sided'"
	//display "`repsnum'"
	//display "`alpha'"

	local state = c(seed)
	set seed `state'

	///////////////////////////////////////////// method monte carlo //////////////////////////////////////////////
	if (strmatch(strlower("`method'"),"mc") )        mata: dbelmc(datmatx,`repsnum',"`varlist'","`sided'")


	///////////////////////////////////////////////// method table ////////////////////////////////////////////////
	//// one sample
	if (strmatch(strlower("`method'"),"interpolate") & `ncolsflag'==1 )   	{
						mata: ncols=cols(datmatx)
						mata: nrows=nrows(datmatx)
						mata: reg=reg(datmatx)

						sysuse critval1s, clear
						quietly putmata crit1s=(*), replace omitmissing
						mata: ts=dbeltable(reg, nrows, crit1s, datmatx, "`varlist'", "`sided'")

						mata: tableexact(reg, nrows, ncols, ts, crit1s, "`sided'")

						drop _all
						getmata (reg)=reg, replace force
						local regflag=reg

						if (`regflag'==1) {
									drop _all
									getmata (sided alpha n cval reg predict)=crit1s, replace force

									gen invn=1/n
									gen n2=n^2
									gen alphahalf=alpha^0.5
									qui reg cval n invn n2 alpha alphahalf

									mata: betas = st_matrix("e(b)")
									mata: decision=getroot(crit1s, betas, ts, "`sided'")
							}
			}

	//// two sample
	if (strmatch(strlower("`method'"),"interpolate") & `ncolsflag'==2 )   	{
						mata: ncols=cols(datmatx)
						mata: nrows=nrows(datmatx)
						mata: reg=reg(datmatx)

						sysuse critval2s, clear
						quietly putmata crit2s=(*), replace omitmissing
						mata: ts=dbeltable2(reg, nrows, crit2s, datmatx, "`varlist'", "`sided'")

						mata: tableexact2(reg, nrows, ncols, ts, crit2s, "`sided'")

						drop _all
						getmata (reg)=reg, replace force
						local regflag=reg

						if (`regflag'==1) {
									drop _all
									getmata (n m cval alpha sided reg predictn predictm)=crit2s, replace force

									gen invn=1/n
									gen n2=n^2
									gen invm=1/m
									gen m2=m^2
									gen alphahalf=alpha^0.5
									qui reg cval n m invn invm n2 m2 alpha alphahalf

									mata: betas = st_matrix("e(b)")
									mata: decision=getroot2(crit2s, betas, ts, "`sided'")
								}
			}

	//// three sample
	if (strmatch(strlower("`method'"),"interpolate") & `ncolsflag'==3 )   	{
						mata: ncols=cols(datmatx)
						mata: nrows=nrows(datmatx)
						mata: reg=reg(datmatx)

						sysuse critval3s, clear
						quietly putmata crit3s=(*), replace omitmissing

						mata: ts=dbeltable3(reg, nrows, crit3s, datmatx, "`varlist'", "`sided'")

						mata: tableexact3(reg, nrows, ncols, ts, crit3s, "`sided'")

						drop _all
						getmata (reg)=reg, replace force
						local regflag=reg

						if (`regflag'==1) {
									drop _all
									getmata (n m p cval alpha sided reg predictn predictm predictp)=crit3s, replace force

									gen invn=1/n
									gen n2=n^2
									gen invm=1/m
									gen m2=m^2
									gen invp=1/p
									gen p2=p^2
									gen alphahalf=alpha^0.5
									qui reg cval n m p invn invm invp n2 m2 p2 alpha alphahalf

									mata: betas = st_matrix("e(b)")
									mata: decision=getroot3(crit3s, betas, ts, "`sided'")
								}
			}





	////////////////////////////////////////////////// method Bayes ////////////////////////////////////////////
	//// one sample
	if (strmatch(strlower("`method'"),"hybrid") & `ncolsflag'==1)   {
						mata: nrows=nrows(datmatx)
						mata: reg=reg(datmatx)

						sysuse critval1s, clear
						quietly putmata crit1s=(*), replace omitmissing
						mata: ts=dbeltable(reg, nrows, crit1s, datmatx, "`varlist'", "`sided'")

						mata: tableexact(reg, nrows, ncols, ts, crit1s, "`sided'")

						drop _all
						getmata (reg)=reg, replace force
						local regflag=reg

						if (`regflag'==1) {
									drop _all
									getmata (sided alpha n cval reg predict)=crit1s, replace force

									gen invn=1/n
									gen n2=n^2
									qui reg cval n invn n2 alpha

									quietly keep in 1
									gen sig0=e(rmse)
									keep predict sig0

									quietly putmata tempprior=(predict sig0), replace

									mata: betas = st_matrix("e(b)")
									mata: prior_param=getprior(betas, tempprior, `alpha')

									mata: dbelbayes(datmatx, prior_param, "`varlist'", "`sided'",`alpha')
								}
			}

	//// two sample
	if (strmatch(strlower("`method'"),"hybrid") & `ncolsflag'==2)   {
						mata: nrows=nrows(datmatx)
						mata: reg=reg(datmatx)

						sysuse critval2s, clear
						quietly putmata crit2s=(*), replace omitmissing
						mata: ts=dbeltable2(reg, nrows, crit2s, datmatx, "`varlist'", "`sided'")

						mata: tableexact2(reg, nrows, ncols, ts, crit2s, "`sided'")

						drop _all
						getmata (reg)=reg, replace force
						local regflag=reg

						if (`regflag'==1) {
									drop _all
									getmata (n m cval alpha sided reg predictn predictm)=crit2s, replace force

									gen invn=1/n
									gen n2=n^2
									gen invm=1/m
									gen m2=m^2
									qui reg cval n m invn invm n2 m2 alpha

									quietly keep in 1
									gen sig0=e(rmse)
									keep predictn predictm sig0

									quietly putmata tempprior=(predictn predictm sig0), replace

									mata: betas = st_matrix("e(b)")
									mata: prior_param=getprior2(betas, tempprior, `alpha')
									mata: dbelbayes(datmatx, prior_param, "`varlist'", "`sided'",`alpha')
								}
			}


	//// three sample
	if (strmatch(strlower("`method'"),"hybrid") & `ncolsflag'==3)   {
						mata: nrows=nrows(datmatx)
						mata: reg=reg(datmatx)

						sysuse critval3s, clear
						quietly putmata crit3s=(*), replace omitmissing
						mata: ts=dbeltable3(reg, nrows, crit3s, datmatx, "`varlist'", "`sided'")

						mata: tableexact3(reg, nrows, ncols, ts, crit3s, "`sided'")

						drop _all
						getmata (reg)=reg, replace force
						local regflag=reg

						if (`regflag'==1) {
									drop _all
									getmata (n m p cval alpha sided reg predictn predictm predictp)=crit3s, replace force

									gen invn=1/n
									gen n2=n^2
									gen invm=1/m
									gen m2=m^2
									gen invp=1/p
									gen p2=p^2
									qui reg cval n m p invn invm invp n2 m2 p2 alpha

									quietly keep in 1
									gen sig0=e(rmse)
									keep predictn predictm predictp sig0

									quietly putmata tempprior=(predictn predictm predictp sig0), replace

									mata: betas = st_matrix("e(b)")
									mata: prior_param=getprior3(betas, tempprior, `alpha')
									mata: dbelbayes(datmatx, prior_param, "`varlist'", "`sided'",`alpha')
								}
			}


end

///////////////////////////////////////// functions here dowwn /////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////// test statistic  functions //////////////////////////////////////////////

/////////////////////One sided one sample///////////////////////////

mata:
	real scalar function wwsim1s(real scalar xn, real scalar yn)
		{
			real colvector x, y
			real scalar log_vnk

			x = sort(rnormal(xn,1,0,1),1)
			y = sort(-x,1)

			log_vnk=log(elr1sw(x,y))
			//log_vnk

			return(log_vnk)
		}
end

mata:
	real scalar function wwt01s(real colvector x, real colvector y)
		{
			real scalar xn, yn, teststat

			xn=rows(x)
			yn=rows(y)

			x=sort(x,1)
			y=sort(y,1)

			teststat=log(elr1sw(x,y))
			//teststat

			return(teststat)
		}
end

mata:
	real scalar function elr1sw(real colvector q1, real colvector qq)
		{

			real colvector m, rm, mm, L, LL, U, UU, qL, qU, F, Fnum, Fnum3
			real rowvector temp, prodFmat
			real matrix Fmat
			real scalar q1n, qqn, elr, Fnum1, Fnum2

			q1n=rows(q1)
			qqn=rows(qq)
			invqt=1/(q1n)

			m=(round(q1n^0.6)::min(( round(q1n^.9),round(q1n/2) )))

			rm=m#J(q1n,1,1)
			mm=J(rows(m),1,1)#(1::q1n)

			L=mm-rm
			LL=(L:*(L:>0))+(L:<=0)
			U=mm+rm
			UU=(U:*(U:<=q1n))+(U:>q1n)*q1n

			qL=q1[LL]
			qU=q1[UU]

		    F=rowmax((mm_relrank(q1,q1n,qU),mm_relrank(qq,qqn,qU)))-rowmax((mm_relrank(q1,q1n,qL),mm_relrank(qq,qqn,qL)))

			F=editvalue(F,0,invqt)

			Fnum1=2*rm
			Fnum2=1/(2*q1n)
			Fnum3=rm:+1
			Fnum=Fnum1:*(1:-(Fnum2*Fnum3))

			F=Fnum:/(q1n*F)
			Fmat=colshape(F,q1n)
			_transpose(Fmat)

			prodFmat=exp(colsum(log(Fmat)))

			elr=min(prodFmat)

			return(elr)
		}
end


////////////////////// two sided one sample //////////////////////

mata:
	real scalar function wwsim(real scalar xn, real scalar yn)
		{
			real colvector x, y, z
			real scalar log_vnk

			x = sort(rnormal(xn,1,0,1),1)
			y = sort(-x,1)
			z = sort((x\y),1)

			log_vnk=log(elrw(x,z))
			//log_vnk

			return(log_vnk)
		}
end

mata:
	real scalar function wwt0(real colvector x, real colvector y)
		{
			real colvector z
			real scalar xn, yn, teststat

			xn=rows(x)
			yn=rows(y)

			x=sort(x,1)
			y=sort(y,1)
			z=sort((x\y),1)

			teststat=log(elrw(x,z))
			//teststat

			return(teststat)
		}
end

mata:
	real scalar function elrw(real colvector q1, real colvector qq)
		{
			real colvector m, rm, mm, L, LL, U, UU, qL, qU, F, Fnum, Fnum3
			real rowvector temp, prodFmat
			real matrix Fmat
			real scalar q1n, qqn, elr, Fnum1, Fnum2

			q1n=rows(q1)
			qqn=rows(qq)
			invqt=1/(q1n)

			m=(round(q1n^0.6)::min(( round(q1n^.9),round(q1n/2) )))

			rm=m#J(q1n,1,1)
			mm=J(rows(m),1,1)#(1::q1n)

			L=mm-rm
			LL=(L:*(L:>0))+(L:<=0)
			U=mm+rm
			UU=(U:*(U:<=q1n))+(U:>q1n)*q1n

			qL=q1[LL]
			qU=q1[UU]

			F=mm_relrank(qq,rows(qq),qU)-mm_relrank(qq,rows(qq),qL)
			F=editvalue(F,0,1/(qqn))

			Fnum1=2*rm
			Fnum2=1/(2*q1n)
			Fnum3=rm:+1
			Fnum=Fnum1:*(1:-(Fnum2*Fnum3))

			F=Fnum:/(q1n*F)

			Fmat=colshape(F,q1n)
			_transpose(Fmat)

			prodFmat=exp(colsum(log(Fmat)))

			elr=min(prodFmat)

			return(elr)
		}
end

////////////////////// two sided two sample test //////////////////////

mata:
	real scalar function qqsim(real scalar xn, real scalar yn)
		{
			real colvector x, y, z
			real scalar log_vnk

			x = sort(runiform(xn,1),1)
			y = sort(runiform(yn,1),1)
			z = sort((x\y),1)

			log_vnk=log(elr(x,z)*elr(y,z))
			//log_vnk

			return(log_vnk)
		}
end

mata:
	real scalar function qqt0(real colvector x, real colvector y)
		{
			real colvector z
			real scalar xn, yn, teststat

			xn=rows(x)
			yn=rows(y)

			x=sort(x,1)
			y=sort(y,1)
			z=sort((x\y),1)

			teststat=log(elr(x,z)*elr(y,z))
			//teststat

			return(teststat)
		}
end

mata:
	real scalar function elr(real colvector q1, real colvector qq)
		{
			real colvector m, rm, mm, L, LL, U, UU, qL, qU, F
			real rowvector temp, prodFmat
			real matrix Fmat
			real scalar q1n, qqn, elr

			q1n=rows(q1)
			qqn=rows(qq)

			m=(round(q1n^0.6)::min(( round(q1n^.9),round(q1n/2) )))

			rm=m#J(q1n,1,1)
			mm=J(rows(m),1,1)#(1::q1n)

			L=mm-rm
			LL=(L:*(L:>0))+(L:<=0)
			U=mm+rm
			UU=(U:*(U:<=q1n))+(U:>q1n)*q1n

			qL=q1[LL]
			qU=q1[UU]

			F=mm_relrank(qq,rows(qq),qU)-mm_relrank(qq,rows(qq),qL)
			F=editvalue(F,0,1/(qqn))

			F=(2*rm):/(q1n*F)

			Fmat=colshape(F,q1n)
			_transpose(Fmat)

			prodFmat=exp(colsum(log(Fmat)))

			elr=min(prodFmat)

			return(elr)
		}
end

/////////////////////One sided two sample///////////////////////////

mata:
	real scalar function qqsim1s(real scalar xn, real scalar yn)
		{
			real colvector x, y
			real scalar log_vnk

			x = sort(runiform(xn,1),1)
			y = sort(runiform(yn,1),1)

			log_vnk=log(elr1sx(x,y)*elr1sy(y,x))
			//log_vnk

			return(log_vnk)
		}
end

mata:
	real scalar function qqt01s(real colvector x, real colvector y)
		{
			real scalar xn, yn, teststat

			xn=rows(x)
			yn=rows(y)

			x=sort(x,1)
			y=sort(y,1)

			teststat=log(elr1sx(x,y)*elr1sy(y,x))
			//teststat

			return(teststat)
		}
end

mata:
	real scalar function elr1sx(real colvector q1, real colvector qq)
		{

			real colvector m, rm, mm, L, LL, U, UU, qL, qU, F
			real rowvector temp, prodFmat
			real matrix Fmat
			real scalar q1n, qqn, elr

			q1n=rows(q1)
			qqn=rows(qq)
			invqt=1/(q1n+qqn)

			m=(round(q1n^0.6)::min(( round(q1n^.9),round(q1n/2) )))

			rm=m#J(q1n,1,1)
			mm=J(rows(m),1,1)#(1::q1n)

			L=mm-rm
			LL=(L:*(L:>0))+(L:<=0)
			U=mm+rm
			UU=(U:*(U:<=q1n))+(U:>q1n)*q1n

			qL=q1[LL]
			qU=q1[UU]

		    F=rowmax((mm_relrank(q1,q1n,qU),mm_relrank(qq,qqn,qU)))-rowmax((mm_relrank(q1,q1n,qL),mm_relrank(qq,qqn,qL)))

			F=editvalue(F,0,invqt)

			F=(2*rm):/(q1n*F)
			Fmat=colshape(F,q1n)
			_transpose(Fmat)

			prodFmat=exp(colsum(log(Fmat)))

			elr=min(prodFmat)

			return(elr)
		}
end

mata:
	real scalar function elr1sy(real colvector q1, real colvector qq)
		{

			real colvector m, rm, mm, L, LL, U, UU, qL, qU, F
			real rowvector temp, prodFmat
			real matrix Fmat
			real scalar q1n, qqn, elr

			q1n=rows(q1)
			qqn=rows(qq)
			invqt=1/(q1n+qqn)

			m=(round(q1n^0.6)::min(( round(q1n^.9),round(q1n/2) )))

			rm=m#J(q1n,1,1)
			mm=J(rows(m),1,1)#(1::q1n)

			L=mm-rm
			LL=(L:*(L:>0))+(L:<=0)
			U=mm+rm
			UU=(U:*(U:<=q1n))+(U:>q1n)*q1n

			qL=q1[LL]
			qU=q1[UU]

		    F=rowmin((mm_relrank(q1,q1n,qU),mm_relrank(qq,qqn,qU)))-rowmin((mm_relrank(q1,q1n,qL),mm_relrank(qq,qqn,qL)))


			F=editvalue(F,0,invqt)

			F=(2*rm):/(q1n*F)
			Fmat=colshape(F,q1n)
			_transpose(Fmat)

			prodFmat=exp(colsum(log(Fmat)))

			elr=min(prodFmat)

			return(elr)
		}
end


////////////////////// three sample test //////////////////////

mata:
	real scalar function qqsim3s(real scalar xn, real scalar yn, real scalar zn)
		{
			real colvector x, y, z, xyz
			real scalar log_vnk

			x = sort(runiform(xn,1),1)
			y = sort(runiform(yn,1),1)
			z = sort(runiform(zn,1),1)
			xyz = sort((x\y\z),1)

			log_vnk=log(elr3s(x,xyz)*elr3s(y,xyz)*elr3s(z,xyz))
			//log_vnk

			return(log_vnk)
		}
end

mata:
	real scalar function qqt03s(real colvector x, real colvector y, real colvector z)
		{
			real colvector xyz
			real scalar xn, yn, zn, teststat

			xn=rows(x)
			yn=rows(y)
			zn=rows(z)

			x=sort(x,1)
			y=sort(y,1)
			z=sort(z,1)
			xyz=sort((x\y\z),1)

			teststat=log(elr3s(x,xyz)*elr3s(y,xyz)*elr3s(z,xyz))
			//teststat

			return(teststat)
		}
end

mata:
	real scalar function elr3s(real colvector q1, real colvector qq)
		{
			real colvector m, rm, mm, L, LL, U, UU, qL, qU, F
			real rowvector temp, prodFmat
			real matrix Fmat
			real scalar q1n, qqn, elr

			q1n=rows(q1)
			qqn=rows(qq)

			m=(round(q1n^0.6)::min(( round(q1n^.9),round(q1n/2) )))

			rm=m#J(q1n,1,1)
			mm=J(rows(m),1,1)#(1::q1n)

			L=mm-rm
			LL=(L:*(L:>0))+(L:<=0)
			U=mm+rm
			UU=(U:*(U:<=q1n))+(U:>q1n)*q1n

			qL=q1[LL]
			qU=q1[UU]

			F=mm_relrank(qq,rows(qq),qU)-mm_relrank(qq,rows(qq),qL)
			F=editvalue(F,0,1/(qqn))

			F=(2*rm):/(q1n*F)

			Fmat=colshape(F,q1n)
			_transpose(Fmat)

			prodFmat=exp(colsum(log(Fmat)))

			elr=min(prodFmat)
			//elr

			return(elr)
		}
end


//////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////// begin and branch by test //////////////////////////////////////////////


//////////////////////////////////////////////// mc ////////////////////////////////////////////////////
mata:
	void dbelmc(real matrix x, real scalar rnum, string dsnlist, string testside)
	{
		string pvaluec

		if (cols(x)>=1) xx=select(x[.,1],x[.,1]:<.)
		if (cols(x)>=2) yy=select(x[.,2],x[.,2]:<.)
		if (cols(x)>=3) zz=select(x[.,3],x[.,3]:<.)

		mc=J(rnum,1,0)

		//////// two sided one sample ////////
		if (strmatch(testside,"one.sample"))
			{
				nn=1
				yy=-xx
				while(nn<=rnum) {
					  mc[nn,1]=wwsim(rows(xx),rows(yy))
					  ++nn
					}
				mc=sort(mc,1)
				t0=wwt0(xx,yy)
				pvalue=(colsum(mc:>=t0))/rnum

				pvaluec=strofreal(pvalue,"%10.4f")
				if (pvalue<0.0001) pvaluec="<0.0001"
				if (pvalue>0.9999) pvaluec=">0.9999"

				test="DBEL (1-sample, 2-sided)"

				printf("{txt}{space 25}{c |}  T0        P-value\n")
				printf("{hline 25}{c +}{hline 25}\n")
				printf("{txt}%24s {c |} {res}%-10.4f {txt}%6s\n",test, t0, pvaluec)
				printf("{txt}{space 25}\n\n\n")
			}

		//////// one sided one sample ////////
		if ((strmatch(testside,"greater") | strmatch(testside,"less")) & cols(x)==1)
			{
				nn=1
				yy=-xx
				while(nn<=rnum) {
					  mc[nn,1]=wwsim1s(rows(xx),rows(yy))
					  ++nn
					}
				mc=sort(mc,1)
				t0=wwt01s(xx,yy)
				pvalue=(colsum(mc:>=t0))/rnum

				pvaluec=strofreal(pvalue,"%10.4f")
				if (pvalue<0.0001) pvaluec="<0.0001"
				if (pvalue>0.9999) pvaluec=">0.9999"

				if (strmatch(testside,"greater")) 		test="DBEL (1-sample, greater)"
				if (strmatch(testside,"less")) 			test="DBEL (1-sample, less)"

				printf("{txt}{space 25}{c |}  T0        P-value\n")
				printf("{hline 25}{c +}{hline 25}\n")
				printf("{txt}%24s {c |} {res}%-10.4f {txt}%6s\n",test, t0, pvaluec)
				printf("{txt}{space 25}\n\n\n")
			}

		//////// two sided two sample ////////
		if (strmatch(testside,"two.sample"))
			{
				nn=1
				while(nn<=rnum) {
					  mc[nn,1]=qqsim(rows(xx),rows(yy))
					  ++nn
					}
				mc=sort(mc,1)
				t0=qqt0(xx,yy)
				pvalue=(colsum(mc:>=t0))/rnum

				pvaluec=strofreal(pvalue,"%10.4f")
				if (pvalue<0.0001) pvaluec="<0.0001"
				if (pvalue>0.9999) pvaluec=">0.9999"

				test="DBEL (2-sample, 2-sided)"

				printf("{txt}{space 25}{c |}  T0        P-value\n")
				printf("{hline 25}{c +}{hline 25}\n")
				printf("{txt}%24s {c |} {res}%-10.4f {txt}%6s\n",test, t0, pvaluec)
				printf("{txt}{space 25}\n\n\n")
			}


		/////// two sample, one sided  ///////
		if ((strmatch(testside,"greater") | strmatch(testside,"less")) & cols(x)==2)
			{
				nn=1
				while(nn<=rnum) {
					  mc[nn,1]=qqsim1s(rows(xx),rows(yy))
					  ++nn
					}
				mc=sort(mc,1)
				t0=qqt01s(xx,yy)
				pvalue=(colsum(mc:>=t0))/rnum

				pvaluec=strofreal(pvalue,"%10.4f")
				if (pvalue<0.0001) pvaluec="<0.0001"
				if (pvalue>0.9999) pvaluec=">0.9999"

				if (strmatch(testside,"greater")) 		test="DBEL (2-sample, greater)"
				if (strmatch(testside,"less")) 			test="DBEL (2-sample, less)"

				printf("{txt}{space 25}{c |}  T0        P-value\n")
				printf("{hline 25}{c +}{hline 25}\n")
				printf("{txt}%24s {c |} {res}%-10.4f {txt}%6s\n",test, t0, pvaluec)
				printf("{txt}{space 25}\n\n\n")
			}
		///////// three sample /////////

		if (strmatch(testside,"three.sample"))
			{
				nn=1
				while(nn<=rnum) {
					  mc[nn,1]=qqsim3s(rows(xx),rows(yy),rows(zz))
					  ++nn
					}
				mc=sort(mc,1)
				t0=qqt03s(xx,yy,zz)
				pvalue=(colsum(mc:>=t0))/rnum

				pvaluec=strofreal(pvalue,"%10.4f")
				if (pvalue<0.0001) pvaluec="<0.0001"
				if (pvalue>0.9999) pvaluec=">0.9999"

				test="DBEL (3-sample)"

				printf("{txt}{space 25}{c |}  T0        P-value\n")
				printf("{hline 25}{c +}{hline 25}\n")
				printf("{txt}%24s {c |} {res}%-10.4f {txt}%6s\n",test, t0, pvaluec)
				printf("{txt}{space 25}\n\n\n")

			}

	}
end


//////////////////////////////////////////////// table ///////////////////////////////////////////////////

////// one sample //////
mata:
real scalar reg( real matrix x)
	{
		nrows=J(1,3,0)
		reg=J(1,1,0)
		if (cols(x)>=3) nrows[1,3]=colsum((x[.,3]:<.))
		if (cols(x)>=2) nrows[1,2]=colsum((x[.,2]:<.))
		if (cols(x)>=1) nrows[1,1]=colsum((x[.,1]:<.))

		if (nrows[1,1]>30 | nrows[1,2]>30 | nrows[1,3]>30 )  reg=1

		return(reg)
	}
end

mata:
real rowvector nrows( real matrix x)
	{
		nrows=J(1,3,0)
		if (cols(x)>=3) nrows[1,3]=colsum((x[.,3]:<.))
		if (cols(x)>=2) nrows[1,2]=colsum((x[.,2]:<.))
		if (cols(x)>=1) nrows[1,1]=colsum((x[.,1]:<.))

		nrows=nrows[1,1],nrows[1,2],nrows[1,3]

		return(nrows)
	}
end

mata:
real scalar dbeltable(real scalar reg, real rowvector nrows, real matrix crit1s, real matrix x, string dsnlist, string testside)
	{
		if (cols(x)>=1) xx=select(x[.,1],x[.,1]:<.)
		if (cols(x)>=2) yy=select(x[.,2],x[.,2]:<.)
		if (cols(x)>=3) zz=select(x[.,3],x[.,3]:<.)

		if ((strmatch(testside,"greater") | strmatch(testside,"less")) & cols(x)==1)
			{
				crit1s=select(crit1s[.,1::13],crit1s[.,1]:==1)
				yy=-xx
				t0=wwt01s(xx,yy)
				if (reg==1) {
								crit1s=select(crit1s[.,1::13],crit1s[.,3]:>=30)
								if (30<nrows[1,1]  & nrows[1,1]<35) 	crit1s=select(crit1s[.,1::4],crit1s[.,5]:==1)
								if (35<=nrows[1,1] & nrows[1,1]<40)		crit1s=select(crit1s[.,1::4],crit1s[.,6]:==1)
								if (40<=nrows[1,1] & nrows[1,1]<50) 	crit1s=select(crit1s[.,1::4],crit1s[.,7]:==1)
								if (50<=nrows[1,1] & nrows[1,1]<60) 	crit1s=select(crit1s[.,1::4],crit1s[.,8]:==1)
								if (60<=nrows[1,1] & nrows[1,1]<80) 	crit1s=select(crit1s[.,1::4],crit1s[.,9]:==1)
								if (80<=nrows[1,1] & nrows[1,1]<100) 	crit1s=select(crit1s[.,1::4],crit1s[.,10]:==1)
								if (100<=nrows[1,1] & nrows[1,1]<150) 	crit1s=select(crit1s[.,1::4],crit1s[.,11]:==1)
								if (150<=nrows[1,1] & nrows[1,1]<200) 	crit1s=select(crit1s[.,1::4],crit1s[.,12]:==1)
								if (200<=nrows[1,1]) 					crit1s=select(crit1s[.,1::4],crit1s[.,13]:==1)

								addinfo=J(rows(crit1s),2,1)
								addinfo[.,1]=addinfo[.,1]:*reg[1,1]
								addinfo[.,2]=addinfo[.,1]:*nrows[1,1]
								crit1s=(crit1s,addinfo)

							}
				if (reg==0) crit1s=select(crit1s[.,1::4],crit1s[.,3]:==nrows[1,1])

			}

		if (strmatch(testside,"one.sample"))
			{
				crit1s=select(crit1s[.,1::13],crit1s[.,1]:==2)
				yy=-xx
				t0=wwt0(xx,yy)
				if (reg==1) {
								crit1s=select(crit1s[.,1::13],crit1s[.,3]:>=30)
								if (30<nrows[1,1]  & nrows[1,1]<35) 	crit1s=select(crit1s[.,1::4],crit1s[.,5]:==1)
								if (35<=nrows[1,1] & nrows[1,1]<40)		crit1s=select(crit1s[.,1::4],crit1s[.,6]:==1)
								if (40<=nrows[1,1] & nrows[1,1]<50) 	crit1s=select(crit1s[.,1::4],crit1s[.,7]:==1)
								if (50<=nrows[1,1] & nrows[1,1]<60) 	crit1s=select(crit1s[.,1::4],crit1s[.,8]:==1)
								if (60<=nrows[1,1] & nrows[1,1]<80) 	crit1s=select(crit1s[.,1::4],crit1s[.,9]:==1)
								if (80<=nrows[1,1] & nrows[1,1]<100) 	crit1s=select(crit1s[.,1::4],crit1s[.,10]:==1)
								if (100<=nrows[1,1] & nrows[1,1]<150) 	crit1s=select(crit1s[.,1::4],crit1s[.,11]:==1)
								if (150<=nrows[1,1] & nrows[1,1]<200) 	crit1s=select(crit1s[.,1::4],crit1s[.,12]:==1)
								if (200<=nrows[1,1]) 					crit1s=select(crit1s[.,1::4],crit1s[.,13]:==1)

								addinfo=J(rows(crit1s),2,1)
								addinfo[.,1]=addinfo[.,1]:*reg[1,1]
								addinfo[.,2]=addinfo[.,1]:*nrows[1,1]
								crit1s=(crit1s,addinfo)

							}
				if (reg==0) crit1s=select(crit1s[.,1::4],crit1s[.,3]:==nrows[1,1])

			}

		return(t0)
	}

end

mata:

void tableexact(real scalar reg,real rowvector nrows, real scalar ncols, real scalar ts, real matrix crit1s, string testside)
	{
		//////// one sample -- 1s & 2s  ////////
		if (ncols==1 & reg==0)
			{
				Tcrit1s=crit1s'
				//Tcrit1s
				if (ts>Tcrit1s[4,1]) pvalue="pvalue<0.01"
				if (Tcrit1s[4,1]>=ts  & ts>Tcrit1s[4,2]) 	pvalue="0.01<=pvalue<0.02"
				if (Tcrit1s[4,2]>=ts  & ts>Tcrit1s[4,3]) 	pvalue="0.02<=pvalue<0.03"
				if (Tcrit1s[4,3]>=ts  & ts>Tcrit1s[4,4]) 	pvalue="0.03<=pvalue<0.04"
				if (Tcrit1s[4,4]>=ts  & ts>Tcrit1s[4,5]) 	pvalue="0.04<=pvalue<0.05"
				if (Tcrit1s[4,5]>=ts  & ts>Tcrit1s[4,6]) 	pvalue="0.05<=pvalue<0.06"
				if (Tcrit1s[4,6]>=ts  & ts>Tcrit1s[4,7]) 	pvalue="0.06<=pvalue<0.07"
				if (Tcrit1s[4,7]>=ts  & ts>Tcrit1s[4,8]) 	pvalue="0.07<=pvalue<0.08"
				if (Tcrit1s[4,8]>=ts  & ts>Tcrit1s[4,9]) 	pvalue="0.08<=pvalue<0.09"
				if (Tcrit1s[4,9]>=ts  & ts>Tcrit1s[4,10]) 	pvalue="0.09<=pvalue<0.10"
				if (Tcrit1s[4,10]>=ts & ts>Tcrit1s[4,11]) 	pvalue="0.10<=pvalue<0.11"
				if (Tcrit1s[4,11]>=ts  & ts>Tcrit1s[4,12]) 	pvalue="0.11<=pvalue<0.12"
				if (Tcrit1s[4,12]>=ts  & ts>Tcrit1s[4,13]) 	pvalue="0.12<=pvalue<0.13"
				if (Tcrit1s[4,13]>=ts  & ts>Tcrit1s[4,14]) 	pvalue="0.13<=pvalue<0.14"
				if (Tcrit1s[4,14]>=ts  & ts>Tcrit1s[4,15]) 	pvalue="0.14<=pvalue<0.15"
				if (Tcrit1s[4,15]>=ts  & ts>Tcrit1s[4,16]) 	pvalue="0.15<=pvalue<0.50"
				if (Tcrit1s[4,16]>=ts  & ts>Tcrit1s[4,17]) 	pvalue="0.50<=pvalue<0.70"
				if (Tcrit1s[4,17]>=ts  & ts>Tcrit1s[4,18]) 	pvalue="0.70<=pvalue<0.80"
				if (Tcrit1s[4,18]>=ts  & ts>Tcrit1s[4,19]) 	pvalue="0.80<=pvalue<0.90"
				if (Tcrit1s[4,19]>=ts) 						pvalue="pvalue>=0.90"

				if (strmatch(testside,"one.sample")) 	test="DBEL (1-sample, 2-sided)"
				if (strmatch(testside,"greater")) 		test="DBEL (1-sample, greater)"
				if (strmatch(testside,"less")) 			test="DBEL (1-sample, less)"

				printf("{txt}{space 25}{c |}  T0                     P-value\n")
				printf("{hline 25}{c +}{hline 34}\n")
				printf("{txt}%24s {c |} {res}%10.0g {txt}%24s\n",test, ts, pvalue)
				printf("{txt}{space 25}\n\n\n")
			}

	}

end

mata:
	void function getoptim(todo,b,c1,c2,c3,crit,g,H)
	{
		crit=(c1+c2*b+c3*b^0.5)^2
	}
end

mata:
	function invreg(x,c1,c2,c3) return(c1+c2*x+c3*x^0.5)
end

mata:

void tablereg(real scalar ts, string testside, real scalar pvalue)
	{
		string test, pvaluec

		pvaluec=strofreal(pvalue,"%10.4f")
		if (pvalue==111) pvaluec="<0.01"
		if (pvalue==999) pvaluec=">0.9"


		if (strmatch(testside,"one.sample")) 	test="DBEL (1-sample, 2-sided)"
		if (strmatch(testside,"greater")) 		test="DBEL (1-sample, greater)"
		if (strmatch(testside,"less")) 			test="DBEL (1-sample, less)"

		printf("{txt}{space 25}{c |}  T0        P-value\n")
		printf("{hline 25}{c +}{hline 25}\n")
		printf("{txt}%24s {c |} {res}%-10.4f {txt}%6s\n",test, ts, pvaluec)
		printf("{txt}{space 25}\n\n\n")

	}
end

mata:
	void function getroot(real matrix crit1s, real rowvector betas, real scalar ts, string sided)
		{
			real scalar c1, c2, c3, testroot, rc, pvalue, max_cval, min_cval

			pvalue=9.9999

			npred=crit1s[1,6]
			max_cval=colmax(crit1s[.,4])
			min_cval=colmin(crit1s[.,4])

			c1=betas[1,6] + betas[1,1]*npred + betas[1,2]*(1/npred)+ betas[1,3]*(npred^2) - ts
			c2=betas[1,4]
			c3=betas[1,5]

			testroot=invreg(0,c1,c2,c3)*invreg(1,c1,c2,c3)

			if (testroot<0) rc=mm_root(pvalue=.,&invreg(),0,1,0,10000,c1,c2,c3)

			if (testroot>=0 & ts<=max_cval & ts>=min_cval) {
					S=optimize_init()
					optimize_init_argument(S, 1, c1)
					optimize_init_argument(S, 2,  c2)
					optimize_init_argument(S, 3, c3)
					optimize_init_evaluator(S, &getoptim())
					optimize_init_which(S,"min")
					optimize_init_params(S,0.9)
					optimize_init_nmsimplexdeltas(S, (10e-4))
					optimize_init_technique(S, "nm")

					pvalue=round(optimize(S),10e-4)
				}

			if (ts>max_cval)  pvalue=111
			if (ts<min_cval)  pvalue=999

			tablereg(ts, sided, pvalue)
		}
end

////// two sample /////

mata:
	void function getroot2(real matrix crit2s, real rowvector betas, real scalar ts, string sided)
		{
			real scalar c1, c2, c3, testroot, rc, pvalue, max_cval, min_cval

			pvalue=9.9999

			npred=crit2s[1,7]
			mpred=crit2s[2,8]
			max_cval=colmax(crit2s[.,3])
			min_cval=colmin(crit2s[.,3])

			c1=betas[1,9] + betas[1,1]*npred + betas[1,2]*mpred + betas[1,3]*(1/npred) + betas[1,4]*(1/mpred) + betas[1,5]*(npred^2) + betas[1,6]*(mpred^2)  - ts
			c2=betas[1,7]
			c3=betas[1,8]

			rc=mm_root(pvalue=.,&invreg(),0,1,0,10000,c1,c2,c3)
			pvalue=round(pvalue,10e-4)

			if (ts<min_cval | pvalue>0.90)  pvalue=999
			if (ts>max_cval | pvalue<0.01)  pvalue=111

			//ts, min_cval, max_cval, pvalue

			tablereg2(ts, sided, pvalue)
		}
end

mata:
real scalar dbeltable2(real scalar reg, real rowvector nrows, real matrix crit2s, real matrix x, string dsnlist, string testside)
	{
		if (cols(x)>=1) xx=select(x[.,1],x[.,1]:<.)
		if (cols(x)>=2) yy=select(x[.,2],x[.,2]:<.)
		if (cols(x)>=3) zz=select(x[.,3],x[.,3]:<.)

		//////////////////////////////////////////// two sample /////////////////////////////////////////////

		if ((strmatch(testside,"greater") | strmatch(testside,"less")) & cols(x)==2)
			{
				crit2s=select(crit2s[.,1::14],crit2s[.,5]:==1) //selects sided
				t0=qqt01s(xx,yy)

				if (reg==1) {
								crit2s=select(crit2s[.,1::14],crit2s[.,1]:>30 :| crit2s[.,2]:>30)

								if ((30<nrows[1,1]  & nrows[1,1]<35)   |  (30<nrows[1,2]  & nrows[1,2]<35))		crit2s=select(crit2s[.,1::14],crit2s[.,6]:==1)
								if ((35<=nrows[1,1] & nrows[1,1]<40)   |  (35<=nrows[1,2] & nrows[1,2]<40))		crit2s=select(crit2s[.,1::14],crit2s[.,7]:==1)
								if ((40<=nrows[1,1] & nrows[1,1]<50)   |  (40<nrows[1,2]  & nrows[1,2]<50))		crit2s=select(crit2s[.,1::14],crit2s[.,8]:==1)
								if ((50<=nrows[1,1] & nrows[1,1]<60)   |  (50<nrows[1,2]  & nrows[1,2]<60))		crit2s=select(crit2s[.,1::14],crit2s[.,9]:==1)
								if ((60<=nrows[1,1] & nrows[1,1]<80)   |  (60<nrows[1,2]  & nrows[1,2]<80))		crit2s=select(crit2s[.,1::14],crit2s[.,10]:==1)
								if ((80<=nrows[1,1] & nrows[1,1]<100)  |  (80<nrows[1,2]  & nrows[1,2]<100))	crit2s=select(crit2s[.,1::14],crit2s[.,11]:==1)
								if ((100<=nrows[1,1] & nrows[1,1]<150) |  (100<nrows[1,2]  & nrows[1,2]<150))	crit2s=select(crit2s[.,1::14],crit2s[.,12]:==1)
								if ((150<=nrows[1,1] & nrows[1,1]<200) |  (1500<nrows[1,2] & nrows[1,2]<200))	crit2s=select(crit2s[.,1::14],crit2s[.,13]:==1)
								if (200<=nrows[1,1] | 200<=nrows[1,2]) 											crit2s=select(crit2s[.,1::14],crit2s[.,14]:==1)

								if (30>=nrows[1,1])	crit2s=select(crit2s[.,1::14],crit2s[.,1]:==nrows[1,1])
								if (30>=nrows[1,2])	crit2s=select(crit2s[.,1::14],crit2s[.,2]:==nrows[1,2])
								if (30<nrows[1,1])	crit2s=select(crit2s[.,1::14],crit2s[.,1]:>=30)
								if (30<nrows[1,2])	crit2s=select(crit2s[.,1::14],crit2s[.,2]:>=30)

								crit2s=crit2s[.,1::5]

								addinfo=J(rows(crit2s),3,1)
								addinfo[.,1]=addinfo[.,1]:*reg[1,1]
								addinfo[.,2]=addinfo[.,1]:*nrows[1,1]
								addinfo[.,3]=addinfo[.,1]:*nrows[1,2]
								crit2s=(crit2s,addinfo)
							}

				if (reg==0) 	crit2s=select(crit2s[.,1::5],crit2s[.,1]:==nrows[1,1] :& crit2s[.,2]:==nrows[1,2])

			}


		if (strmatch(testside,"two.sample") & cols(x)==2)
			{
				crit2s=select(crit2s[.,1::14],crit2s[.,5]:==2)
				t0=qqt0(xx,yy)

				if (reg==1) {
								crit2s=select(crit2s[.,1::14],crit2s[.,1]:>30 :| crit2s[.,2]:>30)

								if ((30<nrows[1,1]  & nrows[1,1]<35)   |  (30<nrows[1,2]  & nrows[1,2]<35))		crit2s=select(crit2s[.,1::14],crit2s[.,6]:==1)
								if ((35<=nrows[1,1] & nrows[1,1]<40)   |  (35<=nrows[1,2] & nrows[1,2]<40))		crit2s=select(crit2s[.,1::14],crit2s[.,7]:==1)
								if ((40<=nrows[1,1] & nrows[1,1]<50)   |  (40<nrows[1,2]  & nrows[1,2]<50))		crit2s=select(crit2s[.,1::14],crit2s[.,8]:==1)
								if ((50<=nrows[1,1] & nrows[1,1]<60)   |  (50<nrows[1,2]  & nrows[1,2]<60))		crit2s=select(crit2s[.,1::14],crit2s[.,9]:==1)
								if ((60<=nrows[1,1] & nrows[1,1]<80)   |  (60<nrows[1,2]  & nrows[1,2]<80))		crit2s=select(crit2s[.,1::14],crit2s[.,10]:==1)
								if ((80<=nrows[1,1] & nrows[1,1]<100)  |  (80<nrows[1,2]  & nrows[1,2]<100))	crit2s=select(crit2s[.,1::14],crit2s[.,11]:==1)
								if ((100<=nrows[1,1] & nrows[1,1]<150) |  (100<nrows[1,2]  & nrows[1,2]<150))	crit2s=select(crit2s[.,1::14],crit2s[.,12]:==1)
								if ((150<=nrows[1,1] & nrows[1,1]<200) |  (1500<nrows[1,2] & nrows[1,2]<200))	crit2s=select(crit2s[.,1::14],crit2s[.,13]:==1)
								if (200<=nrows[1,1] | 200<=nrows[1,2]) 											crit2s=select(crit2s[.,1::14],crit2s[.,14]:==1)

								if (30>=nrows[1,1])	crit2s=select(crit2s[.,1::14],crit2s[.,1]:==nrows[1,1])
								if (30>=nrows[1,2])	crit2s=select(crit2s[.,1::14],crit2s[.,2]:==nrows[1,2])
								if (30<nrows[1,1])	crit2s=select(crit2s[.,1::14],crit2s[.,1]:>=30)
								if (30<nrows[1,2])	crit2s=select(crit2s[.,1::14],crit2s[.,2]:>=30)

								crit2s=crit2s[.,1::5]

								addinfo=J(rows(crit2s),3,1)
								addinfo[.,1]=addinfo[.,1]:*reg[1,1]
								addinfo[.,2]=addinfo[.,1]:*nrows[1,1]
								addinfo[.,3]=addinfo[.,1]:*nrows[1,2]
								crit2s=(crit2s,addinfo)
							}

				if (reg==0) 	crit2s=select(crit2s[.,1::5],crit2s[.,1]:==nrows[1,1] :& crit2s[.,2]:==nrows[1,2])

			}

		return(t0)
	}

end

mata:

void tableexact2(real scalar reg,real rowvector nrows, real scalar ncols, real scalar ts, real matrix crit2s, string testside)
	{
		string test, pvalue

		//////// two sample -- 1s & 2s  ////////
		if (ncols==2 & reg==0)
			{
				Tcrit2s=crit2s'
				if (ts>Tcrit2s[3,1]) pvalue="pvalue<0.01"
				if (Tcrit2s[3,1]>=ts  & ts>Tcrit2s[3,2]) 	pvalue="0.01<=pvalue<0.02"
				if (Tcrit2s[3,2]>=ts  & ts>Tcrit2s[3,3]) 	pvalue="0.02<=pvalue<0.03"
				if (Tcrit2s[3,3]>=ts  & ts>Tcrit2s[3,4]) 	pvalue="0.03<=pvalue<0.04"
				if (Tcrit2s[3,4]>=ts  & ts>Tcrit2s[3,5]) 	pvalue="0.04<=pvalue<0.05"
				if (Tcrit2s[3,5]>=ts  & ts>Tcrit2s[3,6]) 	pvalue="0.05<=pvalue<0.06"
				if (Tcrit2s[3,6]>=ts  & ts>Tcrit2s[3,7]) 	pvalue="0.06<=pvalue<0.07"
				if (Tcrit2s[3,7]>=ts  & ts>Tcrit2s[3,8]) 	pvalue="0.07<=pvalue<0.08"
				if (Tcrit2s[3,8]>=ts  & ts>Tcrit2s[3,9]) 	pvalue="0.08<=pvalue<0.09"
				if (Tcrit2s[3,9]>=ts  & ts>Tcrit2s[3,10]) 	pvalue="0.09<=pvalue<0.10"
				if (Tcrit2s[3,10]>=ts & ts>Tcrit2s[3,11]) 	pvalue="0.10<=pvalue<0.11"
				if (Tcrit2s[3,11]>=ts  & ts>Tcrit2s[3,12]) 	pvalue="0.11<=pvalue<0.12"
				if (Tcrit2s[3,12]>=ts  & ts>Tcrit2s[3,13]) 	pvalue="0.12<=pvalue<0.13"
				if (Tcrit2s[3,13]>=ts  & ts>Tcrit2s[3,14]) 	pvalue="0.13<=pvalue<0.14"
				if (Tcrit2s[3,14]>=ts  & ts>Tcrit2s[3,15]) 	pvalue="0.14<=pvalue<0.15"
				if (Tcrit2s[3,15]>=ts  & ts>Tcrit2s[3,16]) 	pvalue="0.15<=pvalue<0.50"
				if (Tcrit2s[3,16]>=ts  & ts>Tcrit2s[3,17]) 	pvalue="0.50<=pvalue<0.70"
				if (Tcrit2s[3,17]>=ts  & ts>Tcrit2s[3,18]) 	pvalue="0.70<=pvalue<0.80"
				if (Tcrit2s[3,18]>=ts  & ts>Tcrit2s[3,19]) 	pvalue="0.80<=pvalue<0.90"
				if (Tcrit2s[3,19]>=ts) 						pvalue="pvalue>=0.90"

				if (strmatch(testside,"two.sample")) 	test="DBEL (2-sample, 2-sided)"
				if (strmatch(testside,"greater")) 		test="DBEL (2-sample, greater)"
				if (strmatch(testside,"less")) 			test="DBEL (2-sample, less)"

				printf("{txt}{space 25}{c |}  T0                     P-value\n")
				printf("{hline 25}{c +}{hline 34}\n")
				printf("{txt}%24s {c |} {res}%-10.4f {txt}%24s\n",test, ts, pvalue)
				printf("{txt}{space 25}\n\n\n")

			}

	}

end

mata:

void tablereg2(real scalar ts, string testside, real scalar pvalue)
	{
		string test, pvaluec

		pvaluec=strofreal(pvalue,"%10.4f")
		if (pvalue==111) pvaluec="<0.01"
		if (pvalue==999) pvaluec=">0.9"

		if (strmatch(testside,"two.sample")) 	test="DBEL (2-sample, 2-sided)"
		if (strmatch(testside,"greater")) 		test="DBEL (2-sample, greater)"
		if (strmatch(testside,"less")) 			test="DBEL (2-sample, less)"

		printf("{txt}{space 25}{c |}  T0        P-value\n")
		printf("{hline 25}{c +}{hline 25}\n")
		//printf("{txt}%24s {c |} {res}%-10.4f {res}%-10.4f\n",test, ts, pvalue)
		printf("{txt}%24s {c |} {res}%-10.4f {txt}%6s\n",test, ts, pvaluec)
		printf("{txt}{space 25}\n\n\n")

	}
end

////// three sample /////

mata:
	void function getroot3(real matrix crit3s, real rowvector betas, real scalar ts, string sided)
		{
			real scalar c1, c2, c3, testroot, rc, pvalue, max_cval, min_cval

			pvalue=9.9999

			npred=crit3s[1,8]
			mpred=crit3s[1,9]
			ppred=crit3s[1,10]
			max_cval=colmax(crit3s[.,4])
			min_cval=colmin(crit3s[.,4])

			c1=betas[1,12] + betas[1,1]*npred + betas[1,2]*mpred + betas[1,3]*ppred + betas[1,4]*(1/npred) + ///
								betas[1,5]*(1/mpred) + betas[1,6]*(1/ppred) + 	betas[1,7]*(npred^2) + ///
								betas[1,8]*(mpred^2) + betas[1,9]*(ppred^2) - ts
			c2=betas[1,10]
			c3=betas[1,11]

			rc=mm_root(pvalue=.,&invreg(),0,1,0,10000,c1,c2,c3)

			pvalue=round(pvalue,10e-4)

			if (ts<min_cval | pvalue>0.90)  pvalue=999
			if (ts>max_cval | pvalue<0.01)  pvalue=111

			tablereg3(ts, sided, pvalue)
		}
end

mata:
real scalar dbeltable3(real scalar reg, real rowvector nrows, real matrix crit3s, real matrix x, string dsnlist, string testside)
	{
		if (cols(x)>=1) xx=select(x[.,1],x[.,1]:<.)
		if (cols(x)>=2) yy=select(x[.,2],x[.,2]:<.)
		if (cols(x)>=3) zz=select(x[.,3],x[.,3]:<.)

		//////////////////////////////////////////// two sample /////////////////////////////////////////////

		if (strmatch(testside,"three.sample") & cols(x)==3)
			{
				t0=qqt03s(xx,yy,zz)

				if (reg==1) {
								crit3s=select(crit3s[.,1::15],crit3s[.,1]:>30 :| crit3s[.,2]:>30 :|crit3s[.,3]:>30)

								if ((30<nrows[1,1]  & nrows[1,1]<35)   |  (30<nrows[1,2]  & nrows[1,2]<35)	  	|  (30<nrows[1,3]  & nrows[1,3]<35))		crit3s=select(crit3s[.,1::15],crit3s[.,7]:==1)
								if ((35<=nrows[1,1] & nrows[1,1]<40)   |  (35<=nrows[1,2] & nrows[1,2]<40)		|  (35<nrows[1,3]  & nrows[1,3]<40))		crit3s=select(crit3s[.,1::15],crit3s[.,8]:==1)
								if ((40<=nrows[1,1] & nrows[1,1]<50)   |  (40<nrows[1,2]  & nrows[1,2]<50)		|  (40<nrows[1,3]  & nrows[1,3]<50))		crit3s=select(crit3s[.,1::15],crit3s[.,9]:==1)
								if ((50<=nrows[1,1] & nrows[1,1]<60)   |  (50<nrows[1,2]  & nrows[1,2]<60)		|  (50<nrows[1,3]  & nrows[1,3]<60))		crit3s=select(crit3s[.,1::15],crit3s[.,10]:==1)
								if ((60<=nrows[1,1] & nrows[1,1]<80)   |  (60<nrows[1,2]  & nrows[1,2]<80)		|  (60<nrows[1,3]  & nrows[1,3]<80))		crit3s=select(crit3s[.,1::15],crit3s[.,11]:==1)
								if ((80<=nrows[1,1] & nrows[1,1]<100)  |  (80<nrows[1,2]  & nrows[1,2]<100)		|  (80<nrows[1,3]  & nrows[1,3]<100))		crit3s=select(crit3s[.,1::15],crit3s[.,12]:==1)
								if ((100<=nrows[1,1] & nrows[1,1]<150) |  (100<nrows[1,2] & nrows[1,2]<150)		|  (100<nrows[1,3] & nrows[1,3]<150))		crit3s=select(crit3s[.,1::15],crit3s[.,13]:==1)
								if ((150<=nrows[1,1] & nrows[1,1]<200) |  (150<nrows[1,2] & nrows[1,2]<200)		|  (150<nrows[1,3] & nrows[1,3]<200))		crit3s=select(crit3s[.,1::15],crit3s[.,14]:==1)
								if (200<=nrows[1,1] | 200<=nrows[1,2] |  200<=nrows[1,3])																	crit3s=select(crit3s[.,1::15],crit3s[.,15]:==1)

								if (30>=nrows[1,1])	crit3s=select(crit3s[.,1::15],crit3s[.,1]:==nrows[1,1])
								if (30>=nrows[1,2])	crit3s=select(crit3s[.,1::15],crit3s[.,2]:==nrows[1,2])
								if (30>=nrows[1,3])	crit3s=select(crit3s[.,1::15],crit3s[.,3]:==nrows[1,3])
								if (30<nrows[1,1])	crit3s=select(crit3s[.,1::15],crit3s[.,1]:>30)
								if (30<nrows[1,2])	crit3s=select(crit3s[.,1::15],crit3s[.,2]:>30)
								if (30<nrows[1,3])	crit3s=select(crit3s[.,1::15],crit3s[.,3]:>30)

								crit3s=crit3s[.,1::6]

								addinfo=J(rows(crit3s),4,1)
								addinfo[.,1]=addinfo[.,1]:*reg[1,1]
								addinfo[.,2]=addinfo[.,1]:*nrows[1,1]
								addinfo[.,3]=addinfo[.,1]:*nrows[1,2]
								addinfo[.,4]=addinfo[.,1]:*nrows[1,3]
								crit3s=(crit3s,addinfo)

							}

				if (reg==0) 	crit3s=select(crit3s[.,1::6],crit3s[.,1]:==nrows[1,1] :& crit3s[.,2]:==nrows[1,2] :& crit3s[.,3]:==nrows[1,3])

			}
		return(t0)
	}

end

mata:
void tableexact3(real scalar reg,real rowvector nrows, real scalar ncols, real scalar ts, real matrix crit3s, string testside)
	{
		string test, pvalue

		//////// two sample -- 1s & 2s  ////////
		if (ncols==3 & reg==0)
			{
				Tcrit3s=crit3s'
				//Tcrit3s
				//ts
				if (ts>Tcrit3s[4,1]) pvalue="pvalue<0.01"
				if (Tcrit3s[4,1]>=ts  & ts>Tcrit3s[4,2]) 	pvalue="0.01<=pvalue<0.02"
				if (Tcrit3s[4,2]>=ts  & ts>Tcrit3s[4,3]) 	pvalue="0.02<=pvalue<0.03"
				if (Tcrit3s[4,3]>=ts  & ts>Tcrit3s[4,4]) 	pvalue="0.03<=pvalue<0.04"
				if (Tcrit3s[4,4]>=ts  & ts>Tcrit3s[4,5]) 	pvalue="0.04<=pvalue<0.05"
				if (Tcrit3s[4,5]>=ts  & ts>Tcrit3s[4,6]) 	pvalue="0.05<=pvalue<0.06"
				if (Tcrit3s[4,6]>=ts  & ts>Tcrit3s[4,7]) 	pvalue="0.06<=pvalue<0.07"
				if (Tcrit3s[4,7]>=ts  & ts>Tcrit3s[4,8]) 	pvalue="0.07<=pvalue<0.08"
				if (Tcrit3s[4,8]>=ts  & ts>Tcrit3s[4,9]) 	pvalue="0.08<=pvalue<0.09"
				if (Tcrit3s[4,9]>=ts  & ts>Tcrit3s[4,10]) 	pvalue="0.09<=pvalue<0.10"
				if (Tcrit3s[4,10]>=ts & ts>Tcrit3s[4,11]) 	pvalue="0.10<=pvalue<0.11"
				if (Tcrit3s[4,11]>=ts  & ts>Tcrit3s[4,12]) 	pvalue="0.11<=pvalue<0.12"
				if (Tcrit3s[4,12]>=ts  & ts>Tcrit3s[4,13]) 	pvalue="0.12<=pvalue<0.13"
				if (Tcrit3s[4,13]>=ts  & ts>Tcrit3s[4,14]) 	pvalue="0.13<=pvalue<0.14"
				if (Tcrit3s[4,14]>=ts  & ts>Tcrit3s[4,15]) 	pvalue="0.14<=pvalue<0.15"
				if (Tcrit3s[4,15]>=ts  & ts>Tcrit3s[4,16]) 	pvalue="0.15<=pvalue<0.50"
				if (Tcrit3s[4,16]>=ts  & ts>Tcrit3s[4,17]) 	pvalue="0.50<=pvalue<0.70"
				if (Tcrit3s[4,17]>=ts  & ts>Tcrit3s[4,18]) 	pvalue="0.70<=pvalue<0.80"
				if (Tcrit3s[4,18]>=ts  & ts>Tcrit3s[4,19]) 	pvalue="0.80<=pvalue<0.90"
				if (Tcrit3s[4,19]>=ts) 						pvalue="pvalue>=0.90"

				test="DBEL (3-sample)"

				printf("{txt}{space 25}{c |}  T0                     P-value\n")
				printf("{hline 25}{c +}{hline 34}\n")
				printf("{txt}%24s {c |} {res}%-10.4f {txt}%24s\n",test, ts, pvalue)
				printf("{txt}{space 25}\n\n\n")

			}

	}

end

mata:
void tablereg3(real scalar ts, string testside, real scalar pvalue)
	{
		string test, pvaluec

		pvaluec=strofreal(pvalue,"%10.4f")
		if (pvalue==111) pvaluec="<0.01"
		if (pvalue==999) pvaluec=">0.9"


		test="DBEL (3-sample)"

		printf("{txt}{space 25}{c |}  T0        P-value\n")
		printf("{hline 25}{c +}{hline 25}\n")
		printf("{txt}%24s {c |} {res}%-10.4f {txt}%6s\n",test, ts, pvaluec)
		printf("{txt}{space 25}\n\n\n")

	}
end



////////////////////////////////////////////// bayes ////////////////////////////////////////////////

////// one sample /////
mata:
	real matrix function getprior(real matrix betas, real matrix tempprior, real scalar alpha)
		{
			n=tempprior[1,1]
			mu0=betas[1,5] + betas[1,1]*n + betas[1,2]*(1/n) + betas[1,3]*(n^2) ///
						+ betas[1,4]*alpha
			priorparam=(mu0,tempprior[1,2])

			return(priorparam)
		}
end

mata:
	real matrix function getTk(real matrix xx, real scalar tk, string testside)
		{
			mc=J(tk,1,0)
			yy=-xx
			nn=1

			//////// one sample, two sided ////////
			if (strmatch(testside,"one.sample") & cols(xx)==1)
				{
					while(nn<=tk) {
						  mc[nn,1]=wwsim(rows(xx),rows(yy))
						  ++nn
						}
					mc=sort(mc,1)
				}

			//////// one sample, one sided ////////
			if ((strmatch(testside,"greater") | strmatch(testside,"less")) & cols(xx)==1)
				{
					while(nn<=tk) {
						  mc[nn,1]=wwsim1s(rows(xx),rows(yy))
						  ++nn
						}
					mc=sort(mc,1)

				}

			return(mc)
		}
end

mata:
	real matrix function getJk(real scalar tk, real scalar alpha)
		{
			L1=(1-alpha)*tk - (tk^0.5)*ln(tk)
			LL=(1,round(L1))
			U1=(1-alpha)*tk + (tk^0.5)*ln(tk)
			UU=(round(U1),tk)
			Jk=(rowmax(LL), rowmin(UU))

			return(Jk)
		}
end

mata:
	real scalar function getqhatca(real colvector TT, real rowvector JJ, real colvector Jk, real scalar tk, string testside, real scalar u0, real scalar sig0, real scalar alpha)
		{
			Tk=TT[JJ[1,1]::JJ[1,2],1]
			Tkminus=(.\Tk[1::rows(Tk)-1,1])
			//Tk,Tkminus

			exp_part=exp( (-tk/(2*alpha*(1-alpha))):*(((1-alpha):-(Jk:/tk)):^2) )
			integral_pi=normal((Tk:-u0):/sig0) - normal((Tkminus:-u0):/sig0)
			integral_qpi=((sig0/(2*pi()))^0.5):*( exp(-((Tkminus:-u0):^2):/(2*sig0)):-exp(-((Tk:-u0):^2):/(2*sig0))):+(u0:*integral_pi)
			qhatca=colsum(exp_part:*integral_qpi)/colsum(exp_part:*integral_pi)

			return(qhatca)
		}
end

mata:
	real colvector function getf2t2qhat1a(real matrix xx, string testside, real scalar qhat1a)
		{
			t2=200
			T2=getTk(xx, t2, testside)
			sig2=variance(T2)
			h=1.06*(sig2^0.5)*(t2^-0.2)
			f2t2qhat1a=(1/t2)*sqrt(1/(2*pi()*(h^2)))*colsum( exp( (-1/(2*(h^2))):*((qhat1a:-T2):^2)) )

			f2t2qhat1aplus=f2t2qhat1a\T2
			return(f2t2qhat1aplus)
		}
end


mata:
	real matrix function dbelbayes(real matrix x, real rowvector prior_param, string dsnlist, string testside, real scalar alpha)
		{
			if (cols(x)>=1) xx=select(x[.,1],x[.,1]:<.)
			if (cols(x)>=2) yy=select(x[.,2],x[.,2]:<.)
			if (cols(x)>=3) zz=select(x[.,3],x[.,3]:<.)

			if (cols(x)==1) dbelbayes1s(xx, prior_param, testside, alpha)
			if (cols(x)==2) dbelbayes2s(xx, yy,  prior_param, testside, alpha)
			if (cols(x)==3) dbelbayes3s(xx, yy, zz, prior_param, testside, alpha)
		}
end


mata:
	real matrix function dbelbayes1s(real colvector xx, real rowvector prior_param, string testside, real scalar alpha)
		{
			real scalar t0, t1, t2, tt, sig0, u0, V1, qhat1a, f2t2qhat1a, qhatca
			real rowvector JJ
			real colvector T2, J1, TT, f2t2qhat1aplus
			string test, decision

			sig0=prior_param[1,2]^2
			u0=prior_param[1,1]

			t1=200
			JJ=getJk(t1,alpha)
			J1=(JJ[1,1]::JJ[1,2])
			TT=getTk(xx, t1, testside)

			t2=200

			while(t1<=35000) {
						//first cumulative sample
						qhat1a=getqhatca(TT, JJ, J1, t1, testside, u0, sig0, alpha)

						//second sample t2=200
						f2t2qhat1aplus=getf2t2qhat1a(xx, testside, qhat1a)
						f2t2qhat1a=f2t2qhat1aplus[1,1]
						T2=f2t2qhat1aplus[2::rows(f2t2qhat1aplus),1]

						// get variance estimator
						V1=(mm_relrank(T2,rows(T2),qhat1a)*(1-mm_relrank(T2,rows(T2),qhat1a))/((f2t2qhat1a^2)*t2))

						// increment T1
						t1=t1+t2
						tt=t1
						TT=sort(TT\T2,1)
						JJ=getJk(t1,alpha)
						J1=(JJ[1,1]::JJ[1,2])

						if (V1<=sig0) {
								//stop procedure
								t1=999999
							}
						tt
					}

			qhatca=getqhatca(TT, JJ, J1, tt, testside, u0, sig0, alpha)
			yy=-xx
			if (strmatch(testside,"one.sample"))  {
						t0=wwt0(xx,yy)
						test="DBEL (1-sample, 2-sided)"
					}
			if (strmatch(testside,"greater") | strmatch(testside,"less"))	{
						t0=wwt01s(xx,yy)
						if (strmatch(testside,"greater")) 		test="DBEL (1-sample, greater)"
						if (strmatch(testside,"less")) 			test="DBEL (1-sample, less)"
					}

			if (t0>qhatca) 	decision="reject null"
			if (t0<=qhatca) decision="do not reject null"

			printf("{txt}{space 25}{c |}           Posterior Expectation\n")
			printf("{txt}{space 25}{c |}  T0          of Quantiles            Alpha           Decision\n")
			printf("{hline 25}{c +}{hline 70}\n")
			printf("{txt}%24s {c |} {res}%-10.4f %-30.4f %-10.3f {txt}%-24s\n",test, t0, qhatca, alpha, decision)
			printf("{txt}{space 25}\n\n\n")

		}
end

////// two sample //////

mata:
	real matrix function getprior2(real matrix betas, real matrix tempprior, real scalar alpha)
		{
			n=tempprior[1,1]
			m=tempprior[1,2]
			mu0=betas[1,8] + betas[1,1]*n + betas[1,2]*m + betas[1,3]*(1/n) + betas[1,4]*(1/m) + betas[1,5]*(n^2) ///
					+ betas[1,6]*(m^2)+ betas[1,7]*alpha
			priorparam=(mu0,tempprior[1,3])

			return(priorparam)
		}
end

mata:
	real matrix function getTk2(real colvector xx, real colvector yy, real scalar tk, string testside)
		{
			mc=J(tk,1,0)
			nn=1

			//////// two sample, two sided ////////
			if (strmatch(testside,"two.sample"))
				{
					while(nn<=tk) {
						  mc[nn,1]=qqsim(rows(xx),rows(yy))
						  ++nn
						}
					mc=sort(mc,1)
				}

			//////// two sample, one sided ////////
			if (strmatch(testside,"greater") | strmatch(testside,"less"))
				{
					while(nn<=tk) {
						  mc[nn,1]=qqsim1s(rows(xx),rows(yy))
						  ++nn
						}
					mc=sort(mc,1)
				}

			return(mc)
		}
end


mata:
	real colvector function getf2t2qhat1a2(real colvector xx, real colvector yy, string testside, real scalar qhat1a)
		{
			t2=200
			T2=getTk2(xx, yy, t2, testside)
			sig2=variance(T2)
			h=1.06*(sig2^0.5)*(t2^-0.2)
			f2t2qhat1a=(1/t2)*sqrt(1/(2*pi()*(h^2)))*colsum( exp( (-1/(2*(h^2))):*((qhat1a:-T2):^2)) )

			f2t2qhat1aplus=f2t2qhat1a\T2
			return(f2t2qhat1aplus)
		}
end


mata:
	real matrix function dbelbayes2s(real colvector xx, real colvector yy, real rowvector prior_param, string testside, real scalar alpha)
		{
			real scalar t0, t1, t2, tt, sig0, u0, V1, qhat1a, f2t2qhat1a, qhatca
			real rowvector JJ
			real colvector T2, J1, TT, f2t2qhat1aplus
			string test, decision

			sig0=prior_param[1,2]^2
			u0=prior_param[1,1]

			t1=200
			JJ=getJk(t1,alpha)
			J1=(JJ[1,1]::JJ[1,2])
			TT=getTk2(xx, yy, t1, testside)

			t2=200

			while(t1<=35000) {
						//first cumulative sample
						qhat1a=getqhatca(TT, JJ, J1, t1, testside, u0, sig0, alpha)

						//second sample t2=200
						f2t2qhat1aplus=getf2t2qhat1a2(xx, yy, testside, qhat1a)
						f2t2qhat1a=f2t2qhat1aplus[1,1]
						T2=f2t2qhat1aplus[2::rows(f2t2qhat1aplus),1]

						// get variance estimator
						V1=(mm_relrank(T2,rows(T2),qhat1a)*(1-mm_relrank(T2,rows(T2),qhat1a))/((f2t2qhat1a^2)*t2))

						// increment T1
						t1=t1+t2
						tt=t1
						TT=sort(TT\T2,1)
						JJ=getJk(t1,alpha)
						J1=(JJ[1,1]::JJ[1,2])

						if (V1<=sig0) {
								//stop procedure
								t1=999999
							}
						tt
					}

			qhatca=getqhatca(TT, JJ, J1, tt, testside, u0, sig0, alpha)

			if (strmatch(testside,"two.sample"))  {
						t0=qqt0(xx,yy)
						test="DBEL (2-sample, 2-sided)"
					}
			if (strmatch(testside,"greater") | strmatch(testside,"less"))	{
						t0=qqt01s(xx,yy)
						if (strmatch(testside,"greater")) 		test="DBEL (1-sample, greater)"
						if (strmatch(testside,"less")) 			test="DBEL (1-sample, less)"
					}

			if (t0>qhatca) 	decision="reject null"
			if (t0<=qhatca) decision="do not reject null"

			printf("{txt}{space 25}{c |}           Posterior Expectation\n")
			printf("{txt}{space 25}{c |}  T0          of Quantiles            Alpha           Decision\n")
			printf("{hline 25}{c +}{hline 70}\n")
			printf("{txt}%24s {c |} {res}%-10.4f %-30.4f %-10.3f {txt}%-24s\n",test, t0, qhatca, alpha, decision)
			printf("{txt}{space 25}\n\n\n")

		}
end

////// three sample //////

mata:
	real matrix function getprior3(real matrix betas, real matrix tempprior, real scalar alpha)
		{
			n=tempprior[1,1]
			m=tempprior[1,2]
			p=tempprior[1,3]
			mu0=betas[1,11] + betas[1,1]*n + betas[1,2]*m + betas[1,3]*p + betas[1,4]*(1/n) + ///
								betas[1,5]*(1/m) + betas[1,6]*(1/p) + 	betas[1,7]*(n^2) + ///
								betas[1,8]*(m^2) + betas[1,9]*(p^2) + betas[1,10]*alpha
			priorparam=(mu0,tempprior[1,4])

			return(priorparam)
		}
end

mata:
	real matrix function getTk3(real colvector xx, real colvector yy, real colvector zz, real scalar tk, string testside)
		{
			mc=J(tk,1,0)
			nn=1

			//////// three sample ////////
			if (strmatch(testside,"three.sample"))
				{
					while(nn<=tk) {
						  mc[nn,1]=qqsim3s(rows(xx),rows(yy),rows(zz))
						  ++nn
						}
					mc=sort(mc,1)
				}


			return(mc)
		}
end

mata:
	real colvector function getf2t2qhat1a3(real colvector xx, real colvector yy,  real colvector zz, string testside, real scalar qhat1a)
		{
			t2=200
			T2=getTk3(xx, yy, zz, t2, testside)
			sig2=variance(T2)
			h=1.06*(sig2^0.5)*(t2^-0.2)
			f2t2qhat1a=(1/t2)*sqrt(1/(2*pi()*(h^2)))*colsum( exp( (-1/(2*(h^2))):*((qhat1a:-T2):^2)) )

			f2t2qhat1aplus=f2t2qhat1a\T2
			return(f2t2qhat1aplus)
		}
end

mata:
	real matrix function dbelbayes3s(real colvector xx, real colvector yy, real colvector zz, real rowvector prior_param, string testside, real scalar alpha)
		{
			real scalar t0, t1, t2, tt, sig0, u0, V1, qhat1a, f2t2qhat1a, qhatca
			real rowvector JJ
			real colvector T2, J1, TT, f2t2qhat1aplus
			string test, decision

			sig0=prior_param[1,2]^2
			u0=prior_param[1,1]

			t1=200
			JJ=getJk(t1,alpha)
			J1=(JJ[1,1]::JJ[1,2])
			TT=getTk3(xx, yy, zz, t1, testside)

			t2=200

			while(t1<=35000) {
						//first cumulative sample
						qhat1a=getqhatca(TT, JJ, J1, t1, testside, u0, sig0, alpha)

						//second sample t2=200
						f2t2qhat1aplus=getf2t2qhat1a3(xx, yy, zz, testside, qhat1a)
						f2t2qhat1a=f2t2qhat1aplus[1,1]
						T2=f2t2qhat1aplus[2::rows(f2t2qhat1aplus),1]

						// get variance estimator
						V1=(mm_relrank(T2,rows(T2),qhat1a)*(1-mm_relrank(T2,rows(T2),qhat1a))/((f2t2qhat1a^2)*t2))

						// increment T1
						t1=t1+t2
						tt=t1
						TT=sort(TT\T2,1)
						JJ=getJk(t1,alpha)
						J1=(JJ[1,1]::JJ[1,2])

						if (V1<=sig0) {
								//stop procedure
								t1=999999
							}
						tt
					}

			qhatca=getqhatca(TT, JJ, J1, tt, testside, u0, sig0, alpha)

			if (strmatch(testside,"three.sample"))  {
						t0=qqt03s(xx,yy,zz)
						test="DBEL (3-sample)"
					}

			if (t0>qhatca) 	decision="reject null"
			if (t0<=qhatca) decision="do not reject null"

			printf("{txt}{space 25}{c |}           Posterior Expectation\n")
			printf("{txt}{space 25}{c |}  T0          of Quantiles            Alpha           Decision\n")
			printf("{hline 25}{c +}{hline 70}\n")
			printf("{txt}%24s {c |} {res}%-10.4f %-30.4f %-10.3f {txt}%-24s\n",test, t0, qhatca, alpha, decision)
			printf("{txt}{space 25}\n\n\n")


		}
end