DATA_SECTION
  dll_int nobs
  dll_matrix obs(1,nobs,1,2)
 LOC_CALCS
  ofstream ofs("obs");
  ofs << obs  << endl;
PARAMETER_SECTION
  dll_init_bounded_vector pcoff(1,2,.02,1.1);
  dll_init_bounded_vector C1(1,3,-10.0,10.0,2)
  dll_init_bounded_vector C2(1,3,-10.0,10.0,2)
  dll_init_vector mu1(1,2,3)
  dll_init_vector mu2(1,2,3)
  dll_matrix S1(1,2,1,2)
  dll_matrix S2(1,2,1,2)
  dll_vector p(1,2)
  objective_function_value f
PROCEDURE_SECTION
  dvariable psum=sum(pcoff);
  f+=100.*square(log(psum+1.e-20));
  p=pcoff/(psum+1.e-20); // so p’s satisfy constraints
  dvar_matrix tmp1(1,2,1,2);
  dvar_matrix tmp2(1,2,1,2);
  tmp1.initialize();
  tmp2.initialize();
  int ii=1;
  int i=0;
  for (i=1;i<=2;i++) { // fill lower triangle
  for (int j=1;j<=i;j++) {
  tmp1(i,j)=C1(ii);
  tmp2(i,j)=C2(ii);
  ii++;
  }
  }
  S1=tmp1*trans(tmp1); // form S1 S2 from Choleski decomp.
  S2=tmp2*trans(tmp2);
  for (i=1;i<=2;i++) { // to make positive definite
  S1(i,i)+=0.1;
  S2(i,i)+=0.1;
  }
  dvariable det1=sqrt(det(S1));
  dvariable det2=sqrt(det(S2));
  dvar_matrix S1inv=inv(S1);
  dvar_matrix S2inv=inv(S2);
  for (i=1;i<=nobs;i++) // add up minus log-likelihood
  {
  // add the 1.e-10 to avoid log(0) and for robustness
  f-= log(1.e-10+p(1)/det1*exp(-.5*(obs(i)-mu1)*S1inv*(obs(i)-mu1))
  +p(2)/det2*exp(-.5*(obs(i)-mu2)*S2inv*(obs(i)-mu2)));
  }
RUNTIME_SECTION
  maximum_function_evaluations 50,100,10000
REPORT_SECTION
  report << "First mean = " << endl << mu1 << endl;
  report << "First covariance matrix = " << endl<< S1 << endl;
  report << "Second mean = " << endl << mu2 << endl;
  report << "Second covariance matrix = " << endl << S2 << endl;
  
