binormal_ function(x,c=log(10),k1=200,k2=100,tol=1e-4,M=1e6,MM=100,pi1=0.5,pi2=0.5,mu1=-0.5,mu2=0.5,v1=.75,v2=.75)
{

# Test for homogeneity in a two-component normal location/scale mixture. 
# Results are given for the likelihood ratio test, the modified likelihood 
# ratio test (Chen, Chen, and Kalbfleisch), the D-Test based on unpenalized 
# maximum likelihood estimation with and without weighting functions, and the 
# D-Test based on penalized maximum likelihood estimation with and without
# weighting functions.
#
# 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 Input:
# 	x  is a vector of data.
#
# Optional Inputs:
#	c  is a penalty parameter governing the strength of 
#	      a  log(4 p_1 p_2)  term that is attached to
#		the log likelihood function in a penalized maximum
#		likelihood estimation scheme.
#  	k1  is the maximum number of iterations for computing
#		unpenalized maximum likelihood estimates using the 
#		expectation maximization algorithm.
#	k2  is the maximum number of iterations for computing
#		penalized maximum likelihood estimates using the
#		expectation maximization algorithm, with the
#		unpenalized maximum likelihood estimates taken as
#		the starting values.
#	tol  is the largest sum of relative changes in parameter
#		estimates for which the expectation maximization
#		algorithm will terminate short of the maximum
#		number of iterations.
#	M  and  MM  define the parameter space for  mu  as  [-M, M]
#		and for  sigma^2  as  [1/MM, MM].
#     pi1, pi2, mu1, mu2, v1, v2   are initial values for the expectation
#				   maximization algorithm used to obtain
#				   unpenalized maximum likelihood estimates.
#				   These initial values can be obtained
#				   by method of moments estimation or
#				   by visual inspection of an undersmoothed
#				   histogram of the data.
#
# Outputs:
# 	lrt     (Unpenalized) Likelihood Ratio Test:
#		critical values must be estimated via simulation
#	mlrt	Modified (Penalized) Likelihood Ratio Test 
#		developed by Chen, Chen, and Kalbfleisch:
#		critical values must be estimated via simulation
#	dk	D-Test based on unpenalized maximum likelihood estimation:
#		critical values must be estimated via simulation
#     dk.32   with weighting function  exp[-32(x-hat mu_0)^2/hat sigma_0^2]
#     dk.16   with weighting function  exp[-16(x-hat mu_0)^2/hat sigma_0^2]
#     dk.8    with weighting function  exp[- 8(x-hat mu_0)^2/hat sigma_0^2]
#     dkpen   D-Test based on penalized maximum likelihood estimation:
#		critical values must be estimated via simulation
#     dkpen.32   with weighting function  exp[-32(x-hat mu_0)^2/hat sigma_0^2]
#     dkpen.16   with weighting function  exp[-16(x-hat mu_0)^2/hat sigma_0^2]
#     dkpen.8    with weighting function  exp[- 8(x-hat mu_0)^2/hat sigma_0^2]
#	nit1    number of expectation-maximization iterations to obtain
#		the unpenalized maximum likelihood estimates of the 
#		mixture parameters
#	mu1, mu2, pi1, pi2, v1, v2   unpenalized maximum likelihood estimates
#     nit2    number of expectation-maximization iterations to obtain
#		the penalized maximum likelihood estimates of the
#		mixture parameters
#	mu1pen, mu2pen, pi1pen, pi2pen, v1pen, v2pen   
#		penalized maximum likelihood estimates


n_ length(x)

# n  is the "sample size"

mu0_ mean(x)
v0_ (n-1)*var(x)/n

mu1pen_ mu1
mu2pen_ mu2
pi1pen_ pi1
pi2pen_ pi2
v1pen_ v1
v2pen_ v2
lrt_ 0
mlrt_ 0
nit1_ 0
nit2_ 0 

storage.mode(x)_"double"

storage.mode(mu1)_"double"
storage.mode(mu2)_"double"
storage.mode(pi1)_"double"
storage.mode(pi2)_"double"
storage.mode(v1)_"double"
storage.mode(v2)_"double"

storage.mode(mu1pen)_"double"
storage.mode(mu2pen)_"double"
storage.mode(pi1pen)_"double"
storage.mode(pi2pen)_"double"
storage.mode(v1pen)_"double"
storage.mode(v2pen)_"double"

storage.mode(c)_"double"
storage.mode(M)_"double"
storage.mode(MM)_"double"
storage.mode(tol)_"double"

storage.mode(mu0)_"double"
storage.mode(v0)_"double"
storage.mode(lrt)_"double"
storage.mode(mlrt)_"double"


# Here we compute parameter estimates under both
# unpenalized and penalized maximum likelihood schemes.


w_.Fortran("binormal", 
		x=as.double(x),
	        n=as.integer(n), 
	        nit1=as.integer(nit1), 	        
		  nit2=as.integer(nit2), 
		  c=as.double(c),
		  M=as.double(M),
		  MM=as.double(MM),
		  k1=as.integer(k1),
		  k2=as.integer(k2),
		  tol=as.double(tol),
		pi1=as.double(pi1), 
		pi2=as.double(pi2), 
		mu1=as.double(mu1), 
		mu2=as.double(mu2),
		v1=as.double(v1),
		v2=as.double(v2),
		pi1pen=as.double(pi1pen), 
		pi2pen=as.double(pi2pen), 
		mu1pen=as.double(mu1pen), 
		mu2pen=as.double(mu2pen),
		v1pen=as.double(v1pen),
		v2pen=as.double(v2pen),
		mu0=as.double(mu0),
		v0=as.double(v0),
		lrt=as.double(lrt),
		mlrt=as.double(mlrt))

pi1_ w$pi1
pi2_ w$pi2
mu1_ w$mu1
mu2_ w$mu2
v1_ w$v1
v2_ w$v2
lrt_ w$lrt
pi1pen_ w$pi1pen
pi2pen_ w$pi2pen
mu1pen_ w$mu1pen
mu2pen_ w$mu2pen
v1pen_ w$v1pen
v2pen_ w$v2pen
mlrt_ w$mlrt
nit1_ w$nit1
nit2_ w$nit2


# Here we compute the D-test statistics based on 
# the parameter estimates we have.

# unpenalized maximum likelihood

pi0_ -1
sig20_ v0
sig21_ v1
sig22_ v2

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





# penalized maximum likelihood

sig21_ v1pen
sig22_ v2pen

dkp1_ (pi1pen^2 / (2 * sqrt(pi * sig21) ) ) + ( pi2pen^2 / (2 * sqrt(pi * sig22) ) ) + 1 / (2 * sqrt(pi * sig20))
dkp2a_ 2 * pi1pen * exp(-0.5*((mu1pen^2/sig21)+(mu0^2/sig20)) + 0.5*((mu1pen/sig21+mu0/sig20)^2/(sig21^-1+sig20^-1)))
dkp2b_ 2 * pi2pen * exp(-0.5*((mu2pen^2/sig22)+(mu0^2/sig20)) + 0.5*((mu2pen/sig22+mu0/sig20)^2/(sig22^-1+sig20^-1)))
dkp3_ 2 * pi1pen * pi2pen * exp(-0.5*((mu2pen^2/sig22)+(mu1pen^2/sig21)) + 0.5*((mu2pen/sig22+mu1pen/sig21)^2/(sig22^-1+sig21^-1)))
dkpen_ 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 + mu1pen*sig21^-1 - mu0*2*C*sig20^-1
B2.01_ mu0^2*sig20^-1 + mu1pen^2*sig21^-1 - mu0^2*2*C*sig20^-1

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

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

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

B0.22_ sig22^-1 + sig22^-1 - 2*C*sig20^-1
B1.22_ mu2pen*sig22^-1 + mu2pen*sig22^-1 - mu0*2*C*sig20^-1
B2.22_ mu2pen^2*sig22^-1 + mu2pen^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*pi1pen*(1/sqrt(B0.01*2*pi*sig20*sig21))*exp(-0.5*B2.01+0.5*B1.01^2/B0.01)
dktemp02_ pi0*pi2pen*(1/sqrt(B0.02*2*pi*sig20*sig22))*exp(-0.5*B2.02+0.5*B1.02^2/B0.02)
dktemp11_ pi1pen*pi1pen*(1/sqrt(B0.11*2*pi*sig21*sig21))*exp(-0.5*B2.11+0.5*B1.11^2/B0.11)
dktemp12_ pi1pen*pi2pen*(1/sqrt(B0.12*2*pi*sig21*sig22))*exp(-0.5*B2.12+0.5*B1.12^2/B0.12)
dktemp22_ pi2pen*pi2pen*(1/sqrt(B0.22*2*pi*sig22*sig22))*exp(-0.5*B2.22+0.5*B1.22^2/B0.22)

dkpen.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 + mu1pen*sig21^-1 - mu0*2*C*sig20^-1
B2.01_ mu0^2*sig20^-1 + mu1pen^2*sig21^-1 - mu0^2*2*C*sig20^-1

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

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

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

B0.22_ sig22^-1 + sig22^-1 - 2*C*sig20^-1
B1.22_ mu2pen*sig22^-1 + mu2pen*sig22^-1 - mu0*2*C*sig20^-1
B2.22_ mu2pen^2*sig22^-1 + mu2pen^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*pi1pen*(1/sqrt(B0.01*2*pi*sig20*sig21))*exp(-0.5*B2.01+0.5*B1.01^2/B0.01)
dktemp02_ pi0*pi2pen*(1/sqrt(B0.02*2*pi*sig20*sig22))*exp(-0.5*B2.02+0.5*B1.02^2/B0.02)
dktemp11_ pi1pen*pi1pen*(1/sqrt(B0.11*2*pi*sig21*sig21))*exp(-0.5*B2.11+0.5*B1.11^2/B0.11)
dktemp12_ pi1pen*pi2pen*(1/sqrt(B0.12*2*pi*sig21*sig22))*exp(-0.5*B2.12+0.5*B1.12^2/B0.12)
dktemp22_ pi2pen*pi2pen*(1/sqrt(B0.22*2*pi*sig22*sig22))*exp(-0.5*B2.22+0.5*B1.22^2/B0.22)

dkpen.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 + mu1pen*sig21^-1 - mu0*2*C*sig20^-1
B2.01_ mu0^2*sig20^-1 + mu1pen^2*sig21^-1 - mu0^2*2*C*sig20^-1

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

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

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

B0.22_ sig22^-1 + sig22^-1 - 2*C*sig20^-1
B1.22_ mu2pen*sig22^-1 + mu2pen*sig22^-1 - mu0*2*C*sig20^-1
B2.22_ mu2pen^2*sig22^-1 + mu2pen^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*pi1pen*(1/sqrt(B0.01*2*pi*sig20*sig21))*exp(-0.5*B2.01+0.5*B1.01^2/B0.01)
dktemp02_ pi0*pi2pen*(1/sqrt(B0.02*2*pi*sig20*sig22))*exp(-0.5*B2.02+0.5*B1.02^2/B0.02)
dktemp11_ pi1pen*pi1pen*(1/sqrt(B0.11*2*pi*sig21*sig21))*exp(-0.5*B2.11+0.5*B1.11^2/B0.11)
dktemp12_ pi1pen*pi2pen*(1/sqrt(B0.12*2*pi*sig21*sig22))*exp(-0.5*B2.12+0.5*B1.12^2/B0.12)
dktemp22_ pi2pen*pi2pen*(1/sqrt(B0.22*2*pi*sig22*sig22))*exp(-0.5*B2.22+0.5*B1.22^2/B0.22)

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

 
return(lrt,mlrt,dk,dk.32,dk.16,dk.8,dkpen,dkpen.32,dkpen.16,dkpen.8,nit1,mu1,mu2,pi1,pi2,v1,v2,nit2,mu1pen,mu2pen,pi1pen,pi2pen,v1pen,v2pen)
result
}