binormalq_ function(n,pi1,pi2,mu1,mu2,sig21,sig22,mu0,sig20)
{

# Test for homogeneity in a two-component normal location/scale mixture. 
# This is a QUICK test, not requiring the underlying data set itself, 
# but estimates of the mixture parameters must be in hand.
#
# Scenario:
#  	H_0: There exists (mu_0, sigma^2_0) such that
#		p_1 f(x | mu_1, sigma^2_1) + p_2 f(x | mu_2, sigma^2_2)
#	     	may be written as f(x | mu_0, sigma^2_0).
#	H_1: There is no such (mu_0, sigma^2_0).
#
# Required Inputs:
#     n  is the size of the underlying data set ("sample size"). 
#     pi1, pi2, mu1, mu2, sig21, sig22, mu0, sig20  are estimates of the 
#           mixture parameters, which could have been obtained in any way 
#	      (e.g., method of moments, unpenalized maximum likelihood, 
#		penalized maximum likelihood).
#
# Outputs:
#	dk	D-Test statistic, no weighting function
#     dk.32  weighting function  exp[-32(x-hat mu_0)^2/hat sigma_0^2]
#     dk.16  weighting function  exp[-16(x-hat mu_0)^2/hat sigma_0^2]
#     dk.8  weighting function   exp[- 8(x-hat mu_0)^2/hat sigma_0^2]
 

pi0_ -1

dkp1_ (pi1^2 / (2 * sqrt(pi * sig21) ) ) + ( pi2^2 / (2 * sqrt(pi * sig22) ) ) + 1 / (2 * sqrt(pi * sig20))
dkp2a_ 2 * pi1 * exp(-0.5*((mu1^2/sig21)+(mu0^2/sig20)) + 0.5*((mu1/sig21+mu0/sig20)^2/(sig21^-1+sig20^-1)))
dkp2b_ 2 * pi2 * exp(-0.5*((mu2^2/sig22)+(mu0^2/sig20)) + 0.5*((mu2/sig22+mu0/sig20)^2/(sig22^-1+sig20^-1)))
dkp3_ 2 * pi1 * pi2 * exp(-0.5*((mu2^2/sig22)+(mu1^2/sig21)) + 0.5*((mu2/sig22+mu1/sig21)^2/(sig22^-1+sig21^-1)))
dk_ dkp1 - dkp2a/sqrt(2*pi*(sig20+sig21)) - dkp2b/sqrt(2*pi*(sig20+sig22)) + dkp3/sqrt(2*pi*(sig21+sig22))


# weighting function  exp[-32(x - hat mu_0)^2/hat sigma_0^2]

C_ -32

B0.00_ sig20^-1 + sig20^-1 - 2*C*sig20^-1
B1.00_ mu0*sig20^-1 + mu0*sig20^-1 - mu0*2*C*sig20^-1
B2.00_ mu0^2*sig20^-1 + mu0^2*sig20^-1 - mu0^2*2*C*sig20^-1

B0.01_ sig20^-1 + sig21^-1 - 2*C*sig20^-1
B1.01_ mu0*sig20^-1 + mu1*sig21^-1 - mu0*2*C*sig20^-1
B2.01_ mu0^2*sig20^-1 + mu1^2*sig21^-1 - mu0^2*2*C*sig20^-1

B0.02_ sig20^-1 + sig22^-1 - 2*C*sig20^-1
B1.02_ mu0*sig20^-1 + mu2*sig22^-1 - mu0*2*C*sig20^-1
B2.02_ mu0^2*sig20^-1 + mu2^2*sig22^-1 - mu0^2*2*C*sig20^-1

B0.11_ sig21^-1 + sig21^-1 - 2*C*sig20^-1
B1.11_ mu1*sig21^-1 + mu1*sig21^-1 - mu0*2*C*sig20^-1
B2.11_ mu1^2*sig21^-1 + mu1^2*sig21^-1 - mu0^2*2*C*sig20^-1

B0.12_ sig21^-1 + sig22^-1 - 2*C*sig20^-1
B1.12_ mu1*sig21^-1 + mu2*sig22^-1 - mu0*2*C*sig20^-1
B2.12_ mu1^2*sig21^-1 + mu2^2*sig22^-1 - mu0^2*2*C*sig20^-1

B0.22_ sig22^-1 + sig22^-1 - 2*C*sig20^-1
B1.22_ mu2*sig22^-1 + mu2*sig22^-1 - mu0*2*C*sig20^-1
B2.22_ mu2^2*sig22^-1 + mu2^2*sig22^-1 - mu0^2*2*C*sig20^-1



dktemp00_ pi0*pi0*(1/sqrt(B0.00*2*pi*sig20*sig20))*exp(-0.5*B2.00+0.5*B1.00^2/B0.00)
dktemp01_ pi0*pi1*(1/sqrt(B0.01*2*pi*sig20*sig21))*exp(-0.5*B2.01+0.5*B1.01^2/B0.01)
dktemp02_ pi0*pi2*(1/sqrt(B0.02*2*pi*sig20*sig22))*exp(-0.5*B2.02+0.5*B1.02^2/B0.02)
dktemp11_ pi1*pi1*(1/sqrt(B0.11*2*pi*sig21*sig21))*exp(-0.5*B2.11+0.5*B1.11^2/B0.11)
dktemp12_ pi1*pi2*(1/sqrt(B0.12*2*pi*sig21*sig22))*exp(-0.5*B2.12+0.5*B1.12^2/B0.12)
dktemp22_ pi2*pi2*(1/sqrt(B0.22*2*pi*sig22*sig22))*exp(-0.5*B2.22+0.5*B1.22^2/B0.22)

dk.32_ dktemp00+2*dktemp01+2*dktemp02+dktemp11+2*dktemp12+dktemp22


# weighting function  exp[-16(x - hat mu_0)^2/hat sigma_0^2]


C_ -16

B0.00_ sig20^-1 + sig20^-1 - 2*C*sig20^-1
B1.00_ mu0*sig20^-1 + mu0*sig20^-1 - mu0*2*C*sig20^-1
B2.00_ mu0^2*sig20^-1 + mu0^2*sig20^-1 - mu0^2*2*C*sig20^-1

B0.01_ sig20^-1 + sig21^-1 - 2*C*sig20^-1
B1.01_ mu0*sig20^-1 + mu1*sig21^-1 - mu0*2*C*sig20^-1
B2.01_ mu0^2*sig20^-1 + mu1^2*sig21^-1 - mu0^2*2*C*sig20^-1

B0.02_ sig20^-1 + sig22^-1 - 2*C*sig20^-1
B1.02_ mu0*sig20^-1 + mu2*sig22^-1 - mu0*2*C*sig20^-1
B2.02_ mu0^2*sig20^-1 + mu2^2*sig22^-1 - mu0^2*2*C*sig20^-1

B0.11_ sig21^-1 + sig21^-1 - 2*C*sig20^-1
B1.11_ mu1*sig21^-1 + mu1*sig21^-1 - mu0*2*C*sig20^-1
B2.11_ mu1^2*sig21^-1 + mu1^2*sig21^-1 - mu0^2*2*C*sig20^-1

B0.12_ sig21^-1 + sig22^-1 - 2*C*sig20^-1
B1.12_ mu1*sig21^-1 + mu2*sig22^-1 - mu0*2*C*sig20^-1
B2.12_ mu1^2*sig21^-1 + mu2^2*sig22^-1 - mu0^2*2*C*sig20^-1

B0.22_ sig22^-1 + sig22^-1 - 2*C*sig20^-1
B1.22_ mu2*sig22^-1 + mu2*sig22^-1 - mu0*2*C*sig20^-1
B2.22_ mu2^2*sig22^-1 + mu2^2*sig22^-1 - mu0^2*2*C*sig20^-1



dktemp00_ pi0*pi0*(1/sqrt(B0.00*2*pi*sig20*sig20))*exp(-0.5*B2.00+0.5*B1.00^2/B0.00)
dktemp01_ pi0*pi1*(1/sqrt(B0.01*2*pi*sig20*sig21))*exp(-0.5*B2.01+0.5*B1.01^2/B0.01)
dktemp02_ pi0*pi2*(1/sqrt(B0.02*2*pi*sig20*sig22))*exp(-0.5*B2.02+0.5*B1.02^2/B0.02)
dktemp11_ pi1*pi1*(1/sqrt(B0.11*2*pi*sig21*sig21))*exp(-0.5*B2.11+0.5*B1.11^2/B0.11)
dktemp12_ pi1*pi2*(1/sqrt(B0.12*2*pi*sig21*sig22))*exp(-0.5*B2.12+0.5*B1.12^2/B0.12)
dktemp22_ pi2*pi2*(1/sqrt(B0.22*2*pi*sig22*sig22))*exp(-0.5*B2.22+0.5*B1.22^2/B0.22)

dk.16_ dktemp00+2*dktemp01+2*dktemp02+dktemp11+2*dktemp12+dktemp22



# weighting function  exp[-8(x - hat mu_0)^2/hat sigma_0^2]


C_ -8

B0.00_ sig20^-1 + sig20^-1 - 2*C*sig20^-1
B1.00_ mu0*sig20^-1 + mu0*sig20^-1 - mu0*2*C*sig20^-1
B2.00_ mu0^2*sig20^-1 + mu0^2*sig20^-1 - mu0^2*2*C*sig20^-1

B0.01_ sig20^-1 + sig21^-1 - 2*C*sig20^-1
B1.01_ mu0*sig20^-1 + mu1*sig21^-1 - mu0*2*C*sig20^-1
B2.01_ mu0^2*sig20^-1 + mu1^2*sig21^-1 - mu0^2*2*C*sig20^-1

B0.02_ sig20^-1 + sig22^-1 - 2*C*sig20^-1
B1.02_ mu0*sig20^-1 + mu2*sig22^-1 - mu0*2*C*sig20^-1
B2.02_ mu0^2*sig20^-1 + mu2^2*sig22^-1 - mu0^2*2*C*sig20^-1

B0.11_ sig21^-1 + sig21^-1 - 2*C*sig20^-1
B1.11_ mu1*sig21^-1 + mu1*sig21^-1 - mu0*2*C*sig20^-1
B2.11_ mu1^2*sig21^-1 + mu1^2*sig21^-1 - mu0^2*2*C*sig20^-1

B0.12_ sig21^-1 + sig22^-1 - 2*C*sig20^-1
B1.12_ mu1*sig21^-1 + mu2*sig22^-1 - mu0*2*C*sig20^-1
B2.12_ mu1^2*sig21^-1 + mu2^2*sig22^-1 - mu0^2*2*C*sig20^-1

B0.22_ sig22^-1 + sig22^-1 - 2*C*sig20^-1
B1.22_ mu2*sig22^-1 + mu2*sig22^-1 - mu0*2*C*sig20^-1
B2.22_ mu2^2*sig22^-1 + mu2^2*sig22^-1 - mu0^2*2*C*sig20^-1



dktemp00_ pi0*pi0*(1/sqrt(B0.00*2*pi*sig20*sig20))*exp(-0.5*B2.00+0.5*B1.00^2/B0.00)
dktemp01_ pi0*pi1*(1/sqrt(B0.01*2*pi*sig20*sig21))*exp(-0.5*B2.01+0.5*B1.01^2/B0.01)
dktemp02_ pi0*pi2*(1/sqrt(B0.02*2*pi*sig20*sig22))*exp(-0.5*B2.02+0.5*B1.02^2/B0.02)
dktemp11_ pi1*pi1*(1/sqrt(B0.11*2*pi*sig21*sig21))*exp(-0.5*B2.11+0.5*B1.11^2/B0.11)
dktemp12_ pi1*pi2*(1/sqrt(B0.12*2*pi*sig21*sig22))*exp(-0.5*B2.12+0.5*B1.12^2/B0.12)
dktemp22_ pi2*pi2*(1/sqrt(B0.22*2*pi*sig22*sig22))*exp(-0.5*B2.22+0.5*B1.22^2/B0.22)

dk.8_ dktemp00+2*dktemp01+2*dktemp02+dktemp11+2*dktemp12+dktemp22


return(dk,dk.32,dk.16,dk.8)
result
}