compute M = 50
    compute N = 5
    all M
    dec vect v(M)  simul1(M)  simul2(M)  u1(M)  u2(M)  v1(M)  v2(M)  j1(M)  j2(M)  j3(M) 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
    
    compute muxs = 0.049969029, muTs = 0.0, sigmax1s = 1.439360726, sigmaT1s = 1.0, sigmaTx1s = 0.0, sigmax2s = 1.439360726, sigmaT2s = 1.0, sigmaTx2s = 0.0

    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))
    nonlin mux muT sigmax1 sigmaT1 sigmaTx1 sigmax2 sigmaT2 sigmaTx2
    compute k=0
     while k<=N {
       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)

       frml y = $
             g1  =(1/sigmax1)*%density((x(t)-mux)/sigmax1),$
             g21 = -((1-((sigmax1)^(-2)*sigmaTx1))*(x(t)-mux))/sqrt(sigmaT1^2-(sigmax1^(-2)*sigmaTx1^2)),$
             g22 = (muT-mux)/sqrt(sigmaT1^2-(sigmax1^(-2)*(sigmaTx1)^2)),$
             g2  = %cdf(g21+g22),$
             g3  = %cdf((muT-mux)/sqrt(sigmaT1^2+(sigmax1)^2+(-2)*sigmaTx1)),$
             g4  = %cdf((mux-muT)/sqrt(sigmaT2^2+(sigmax2)^2 +(-2)*sigmaTx2)),$
             h1  = (1/sigmax2)*%density((x(t)-mux)/sigmax2),$
             h21 = ((1-((sigmax2)^(-2)*sigmaTx2))*(x(t)-mux))/sqrt((sigmaT2^2-(sigmax2^(-2)*sigmaTx2^2))),$
             h22 = (mux-muT)/sqrt(sigmaT2^2-(sigmax2^(-2)*sigmaTx2^2)),$
             h2  = %cdf(h21+h22),$
             log((g1*g2)+(h1*h2)) - log(g3+g4)
             compute mux = 0.0 ,  muT = 0.05 ,  sigmax1 = 1.5 , sigmaT1 =0.95 ,  sigmaTx1 = 0.01 , sigmax2 = 1.5, sigmaT2 = 0.99  , sigmaTx2 = 0.01
             maximize(pmethod=simplex,piters = 9,subiters =1000, iters = 5000, method=bfgs, cvcrit=0.000001, noprint) y
             
              report(action = define, hlabels= || 'mux', 'muT', 'sigmax1', 'sigmaT1', 'sigmaTx1', 'sigmax2', 'sigmaT2', 'sigmaTx2'||)
              if %converged==1 {
              compute k=k+1
              report(row = new, atcol = 1) mux muT sigmax1 sigmaT1 sigmaTx1 sigmax2 sigmaT2 sigmaTx2
              com mux_hat(k) = mux,  muT_hat(k) = muT
              com sigmax1_hat(k) = sigmax1, sigmaT1_hat(k) = sigmaT1
              com sigmax2_hat(k) = sigmaTx1, sigmaT2_hat(k) = sigmax2
              com sigmaTx1_hat(k) = sigmaT2, sigmaTx2_hat(k) = sigmaTx2
               

              report(row=current, atcol=1, align=decimal) mux_hat(k)
              report(row=current, atcol=2, align = decimal) muT_hat(k)
              report(row = current, atcol=3, align = decimal) sigmax1_hat(k)
              report(row = current, atcol=4,  align = decimal) sigmaT1_hat(k)
              report(row= current, atcol=5,  align = decimal) sigmaTx1_hat(k)
              report(row = current, atcol=6, align = decimal) sigmax2_hat(k)
              report(row = current, atcol=7,  align = decimal) sigmaT2_hat(k)
              report(row = current, atcol=8,  align = decimal) sigmaTx2_hat(k) 
              report(action = format, align = center, width = 9)
               
              open copy file_50.xls
              report(action=show,window = 'Estimation')
              report(action=show,format=xls,unit=copy)
               } 
               }
                end do while