/* nag_binary_factor (g11sac) Example Program.

 *

 * Copyright 2002 Numerical Algorithms Group.

 *

 * Mark 7, 2002.

 * Mark 7b revised, 2004.

 */



#include <stdio.h>

#include <nag.h>

#include <nag_stdlib.h>

#include <nagg11.h>



int main(void)

{

  /* Scalars */

  double cgetol, chi, rlogl, siglev;

  Integer exit_status, i, pdcm, idf, p, iprint, is,

    j, maxit, n, niter, nrx, pdx;

  NagError fail;

  Nag_OrderType order;

  Nag_Boolean chisqr, gprob;

  char flag;



  /* Arrays */

  double *a = 0, *alpha = 0, *c = 0, *cm = 0, *exf = 0, *expp = 0,

    *g = 0, *obs = 0, *pigam = 0, *xl = 0, *y = 0;

  Integer *iob = 0, *irl = 0;

  Nag_Boolean *x = 0;



#ifdef NAG_COLUMN_MAJOR

#define X(I,J) x[(J-1)*pdx + I - 1]

#define CM(I,J) cm[(J-1)*pdcm + I - 1]

  order = Nag_ColMajor;

#else

#define X(I,J) x[(I-1)*pdx + J - 1]

#define CM(I,J) cm[(I-1)*pdcm + J - 1]

  order = Nag_RowMajor;

#endif



  INIT_FAIL(fail);

  exit_status = 0;

  Vprintf("nag_binary_factor (g11sac) Example Program Results\n");



  /* Skip heading in data file */

  Vscanf("%*[^\n] ");



  Vscanf("%ld%ld%ld%*[^\n] ", &p, &n, &is);

  if (p > 0 && is >= 0)

    {

      /* Allocate arrays */

      pdcm = 2*p;

      nrx = is;

      if ( !(a = NAG_ALLOC(p, double)) ||

           !(alpha = NAG_ALLOC(p, double)) ||

           !(c = NAG_ALLOC(p, double)) ||

           !(cm = NAG_ALLOC(pdcm * 2*p, double)) ||

           !(exf = NAG_ALLOC(is, double)) ||

           !(expp = NAG_ALLOC(p * p, double)) ||

           !(g = NAG_ALLOC(2*p, double)) ||

           !(obs = NAG_ALLOC(p * p, double)) ||

           !(pigam = NAG_ALLOC(p, double)) ||

           !(xl = NAG_ALLOC(is, double)) ||

           !(y = NAG_ALLOC(is, double)) ||

           !(iob = NAG_ALLOC(is, Integer)) ||

           !(irl = NAG_ALLOC(is, Integer)) ||

           !(x = NAG_ALLOC(nrx * p, Nag_Boolean)) )

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }



      if (order == Nag_ColMajor)

        pdx = nrx;

      else

        pdx = p;



      for (i = 1; i <= is; ++i)

        {

          Vscanf("%ld", &irl[i-1]);

          for (j = 1; j <= p; ++j)

            {

              Vscanf(" %c", &flag);

              X(i,j) = (flag == 'T');

            }

          Vscanf("%*[^\n] ");

        }

      gprob = Nag_FALSE;

      for (i = 1; i <= p; ++i)

        {

          a[i-1] = 0.5;

          c[i-1] = 0.0;

        }



      /* Set iprint > 0 to obtain intermediate output */

      iprint = -1;

      cgetol = 1e-4;

      maxit = 1000;

      chisqr = Nag_TRUE;



      /* nag_binary_factor (g11sac).

       * Contingency table, latent variable model for binary data

       */

      nag_binary_factor(order, p, n, gprob, is, x, pdx, irl, a, c, iprint, 0,

                        cgetol, maxit, chisqr, &niter, alpha, pigam, cm,

                        pdcm, g, expp, p, obs, exf, y, iob, &rlogl, &chi,

                        &idf, &siglev,  &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_binary_factor (g11sac).\n%s\n", fail.message);

          exit_status = 1;

          goto END;

        }



      Vprintf("\n");

      Vprintf("Item     Alpha      (s.e.)          Pi       (s.e.)\n");

      for (i=1; i<=p; i++)

        Vprintf("  %ld     %g  (%10g)   %g  (%10g)\n", i, alpha[i-1],

                CM(2*i-1,2*i-1),

                pigam[i-1], CM(2*i,2*i));

      Vprintf("\n");

      Vprintf("Index   Observed   Expected     Theta    Pattern\n");

      Vprintf("        Frequency  Frequency    Score\n");

      for (i=1; i<=is; i++)

        {

          Vprintf("  %2ld       %3ld      %7g  %10g    ", i,

                  irl[i-1], exf[i-1], y[i-1]);

          for (j=1; j<=p; j++)

            Vprintf("%s",X(i,j)==1?"T":"F");

          Vprintf("\n");

        }

      Vprintf("\n");

      Vprintf("Chi-squared test statistic = %g\n", chi);

      Vprintf("Degrees of freedom =         %ld\n",idf);

      Vprintf("Significance =               %g\n",siglev);

    }



 END:

  if (a) NAG_FREE(a);

  if (alpha) NAG_FREE(alpha);

  if (c) NAG_FREE(c);

  if (cm) NAG_FREE(cm);

  if (exf) NAG_FREE(exf);

  if (expp) NAG_FREE(expp);

  if (g) NAG_FREE(g);

  if (obs) NAG_FREE(obs);

  if (pigam) NAG_FREE(pigam);

  if (xl) NAG_FREE(xl);

  if (y) NAG_FREE(y);

  if (iob) NAG_FREE(iob);

  if (irl) NAG_FREE(irl);

  if (x) NAG_FREE(x);



  return exit_status;

}