compute M = 50
compute N = 1000
all M
dec vect v(M)  simul1(M)  simul2(M)  u1(M)  u2(M)  v1(M)  v2(M)  j1(M)  j2(M)  j3(N) s(M)  x(M)
dec real  mux  muT  sigmax1  sigmaT1  sigmaTx1  sigmax2   sigmaT2  sigmaTx2
set mux_hat 1 N = 0; set muT_hat 1 N = 0; set sigmax1_hat 1 N = 0; set sigmaT1_hat 1 N = 0; set sigmaTx1_hat 1 N = 0
set sigmax2_hat 1 N = 0; set sigmaT2_hat 1 N = 0;  set sigmaTx2_hat 1 N = 0
set bias_1 1 N = 0; set bias_2  1  N = 0.; set bias_3 1 N = 0; set bias_4 1 N = 0; set bias_5 1 N = 0
set bias_6 1 N = 0; set bias_7  1 N = 0.;  set bias_8 1 N = 0
set mse_1 1 N = 0; set mse_2  1  N = 0; set mse_3 1 N = 0; set mse_4 1 N = 0; set mse_5 1 N = 0
set mse_6 1 N = 0; set mse_7  1 N = 0;  set mse_8 1 N = 0
compute muxs = 0.0, muTs = 1.0, sigmax1s = 1.0, sigmaT1s = 0.5, sigmaTx1s = 0.45, sigmax2s = 1.0, sigmaT2s = 1.0, sigmaTx2s = 0.5

do k = 1,N
compute lambda01s = (muTs-muxs)/sqrt(sigmaT1s^2-((sigmax1s)^(-2)*sigmaTx1s^2))
compute lambda02s = (muTs-muxs)/sqrt(sigmaT2s^2-((sigmax2s)^(-2)*sigmaTx2s^2))
compute lambda11s = -sigmax1s*(1-(sigmax1s)^(-2)*sigmaTx1s)/sqrt(sigmaT1s^2-((sigmax1s)^(-2)*(sigmaTx1s)^2))
compute lambda12s = -sigmax2s*(1-(sigmax2s)^(-2)*sigmaTx2s)/sqrt(sigmaT2s^2-((sigmax2s)^(-2)*(sigmaTx2s)^2))
compute as = -lambda01s/sqrt(1+lambda11s^2)
compute bs = -lambda02s/sqrt(1+lambda12s^2)
compute cs = 1/sqrt(1+lambda11s^2)
compute ds = 1/sqrt(1+lambda12s^2)
compute g3s = (muTs-muxs)/sqrt(sigmaT1s^2+sigmax1s^2+(-2)*sigmaTx1s)
compute g4s = (muxs-muTs)/sqrt(sigmaT2s^2+sigmax2s^2+(-2)*sigmaTx2s)
compute ps = %cdf(g3s)/(%cdf(g3s)+%cdf(g4s))
ewise v(i)= %ran(1.0)
ewise simul1(i)= %uniform(0,1)
ewise u1(i) = %invnormal(%cdf(as)+((1-%cdf(as))*simul1(i)))
ewise u2(i) = %invnormal(%cdf(bs)*simul1(i))
ewise simul2(i)= %uniform(0,1)
ewise j1(i) = %sign(simul2(i)-ps)
ewise j2(i) = j1(i)+1
ewise j3(i) = %sign(j2(i))
ewise s(i) = cs*(lambda11s*u1(i)+v(i))*(1-j3(i)) + ds*(lambda12s*u2(i)+v(i))*j3(i)
ewise x(i) = (sigmax1s*s(i)+muxs)*(1-j3(i)) + (sigmax2s*s(i)+muxs)*j3(i)

do i=1,M
nonlin mux muT sigmax1 sigmaT1 sigmaTx1 sigmax2 sigmaT2  sigmaTx2
frml g1 = (1/sigmax1)*%density((x(i)-mux)/sigmax1)
frml g21 = -((1-((sigmax1)^(-2)*sigmaTx1))*(x(i)-mux))/sqrt(sigmaT1^2-(sigmax1^(-2)*sigmaTx1^2))
frml g22 = (muT-mux)/sqrt(sigmaT1^2-(sigmax1^(-2)*(sigmaTx1)^2))
frml g2 = %cdf(g21+g22)
frml g3 = %cdf((muT-mux)/sqrt(sigmaT1^2+(sigmax1)^2+(-2)*sigmaTx1))
frml g4 = %cdf((mux-muT)/sqrt(sigmaT2^2+(sigmax2)^2 +(-2)*sigmaTx2))
frml h1 = (1/sigmax2)*%density((x(i)-mux)/sigmax2)
frml h21 = ((1-((sigmax2)^(-2)*sigmaTx2))*(x(i)-mux))/sqrt((sigmaT2^2-(sigmax2^(-2)*sigmaTx2^2)))
frml h22 = (mux-muT)/sqrt(sigmaT2^2-(sigmax2^(-2)*sigmaTx2^2))
frml h2 = %cdf(h21+h22)
frml y = log((g1*g2)+(h1*h2)) - log(g3+g4)
end do i

compute mux = 0.0,  muT = 1.0,  sigmax1 = 1.0, sigmaT1 = 0.5,  sigmaTx1 = 0.45, sigmax2 = 1.0, sigmaT2 = 1.0, sigmaTx2 = 0.5
maximize(method=simplex,iters = 9) y
maximize(iters = 5000, method=bfgs, cvcrit=0.0001) y

com mux_hat(k) = %beta(1),  muT_hat(k) = %beta(2)
com sigmax1_hat(k) = %beta(3), sigmaT1_hat(k) = %beta(4)
com sigmax2_hat(k) = %beta(5), sigmaT2_hat(k) = %beta(6)
com sigmaTx1_hat(k) = %beta(7), sigmaTx2_hat(k) = %beta(8)

com bias_1(k) = muxs - mux_hat(k),  bias_2(k) = muTs - muT_hat(k)
com bias_3(k) = sigmax1s - sigmax1_hat(k),  bias_4(k) = sigmaT1s - sigmaT1_hat(k)
com bias_5(k) = sigmax2s - sigmax2_hat(k), bias_6(k) = sigmaT2s - sigmaT2_hat(k)
com bias_7(k) = sigmaTx1s - sigmaTx1_hat(k),  bias_8(k) = sigmaTx2s - sigmaTx2_hat(k)

com mse_1(k) = (bias_1(k))^2, mse_2(k) = (bias_2(k))^2
com mse_3(k) = (bias_3(k))^2, mse_4(k) = (bias_4(k))^2, mse_5(k) = (bias_5(k))^2
com mse_6(k) = (bias_6(k))^2, mse_7(k) = (bias_7(k))^2, mse_8(k) = (bias_8(k))^2

end do k
table / bias_1 bias_2  bias_3  bias_4  bias_5  bias_6  bias_7  bias_8
table / mse_1  mse_2  mse_3  mse_4  mse_5  mse_6  mse_7  mse_8