DATA_SECTION
  int n
 !! n=500;
  matrix x(1,n,1,3)
  number tau
 !! tau=2.0;
  matrix sigma(1,2,1,2)
  matrix ch(1,2,1,2)
  vector beta(1,2)
  matrix mu(1,n,1,3)
  matrix Y(1,n,1,3)
  matrix b(1,n,1,2)
  matrix prob(1,n,1,3)
 LOC_CALCS
  beta=1.0;
  double tau2=square(tau);
  sigma(1,1)=1.0*tau2;
  sigma(2,2)=1.0*tau2;
  sigma(1,2)=0.25*tau2;
  sigma(2,1)=0.25*tau2;
  ch=choleski_decomp(sigma);
  int iseed;
  {
    ifstream ifs("seed");
    ifs >> iseed;
    iseed+=2;
  }
  {
    ofstream ofs("seed");
    ofs << iseed;
  }
  random_number_generator rng(iseed); 
  x.fill_randn(rng);
  b.fill_randn(rng);
  prob.fill_randu(rng);
  for (int i=1;i<=n;i++)
  {
    b(i)=ch*b(i);
    dvector tmp=beta(1)+b(i,1)
      +(beta(2)+b(i,2))*x(i);
    mu(i)=1.0/(1.0+exp(tmp));
    for (int t=1;t<=3;t++)
    {
      if (prob(i,t)<mu(i,t)) 
        Y(i,t)=1.0;
      else
        Y(i,t)=0.0;
    }
  }
  ofstream ofs("analyzer.dat");
  ofs << "# n" << endl;
  ofs << n << endl;
  ofs << "# x" << endl;
  ofs << x << endl;
  ofs << "# Y" << endl;
  ofs << Y << endl;
  cout << Y << endl;
  cout << setprecision(3) << setfixed() <<setw(6) << mu << endl;
  cout << setprecision(3) << setfixed() <<setw(6) << prob << endl;
  exit(1);

PARAMETER_SECTION
  objective_function_value f;
PROCEDURE_SECTION
