mvnorm.c

数据挖掘中的bayes算法,很好的代码

    r1->cv   = cnt    *vals[i]; /* precompute number of cases * value */
    r1->cv2  = r1->cv *vals[i]; /* and number of cases *value^2 */
    r1->cnt += cnt;             /* update the terms that are needed */
    r1->sv  += r1->cv;          /* to compute the expected value */
    r1->sv2 += r1->cv2;         /* and the variance */
    for (e = r1->elems +(k = i); --k >= 0; ) {
      if (vals[k] <= 0)         /* traverse the preceding rows */
        continue;               /* but skip those rows, for */
      r2 = mvn->rows[k];        /* which the value is not positive */
      (--e)->cnt += cnt;
      e->sr  += r1->cv; e->sr2 += r1->cv2;
      e->sc  += r2->cv; e->sc2 += r2->cv2;
      e->src += r1->cv *vals[k];
    }                           /* update the terms needed */
  }                             /* to compute the covariance */
}  /* mvn_addx() */

/*--------------------------------------------------------------------*/

int mvn_calc (MVNORM *mvn, int mode)
{                               /* --- calc. parameters from data */
  int     i, k;                 /* loop variables */
  MVNROW  *row;                 /* to traverse the matrix rows */
  MVNELEM *e;                   /* to traverse the matrix elements */
  double  *x, *v, *c;           /* to traverse the estim. parameters */
  double  cnt, t;               /* number of cases (-1), buffer */
  int     err = 0;              /* return code */

  assert(mvn);                  /* check the function argument */
  if (mode & MVN_EXPVAR) {      /* expected values and variances */
    for (x = mvn->exps +(i = mvn->size); --i >= 0; ) {
      row  = mvn->rows[i];      /* traverse the statistics rows */
      *--x = (row->cnt > 0)     /* compute the expected value */
           ?  row->sv /row->cnt : 0;
      t    = row->sv2 -*x *row->sv;
      cnt  = (mode & MVN_MAXLLH) ? row->cnt : (row->cnt -1);
      mvn->covs[i][i] = ((t > 0) && (cnt > 0)) ? t /cnt : 0;
    }                           /* compute the variance */
  }
  if (mode & MVN_COVAR) {       /* if to compute covariances */
    for (x = mvn->exps, i = 0; i < mvn->size; i++) {
      row = mvn->rows[i];       /* traverse the statistics rows */
      v   = mvn->covs[i] +i;    /* and the covariance rows */
      if (*v <= 0) {            /* if the variance is zero */
        for (k = i; --k >= 0; ) *--v = 0;
        continue;               /* set all covariances to zero */
      }                         /* (to avoid inconsistencies) */
      for (e = row->elems +(k = i); --k >= 0; ) {
        cnt  = (mode & MVN_MAXLLH) ? (--e)->cnt : ((--e)->cnt -1);
        *--v = (cnt > 0)        /* compute the covariance */
             ? (e->src -x[k] *e->sr -x[i] *e->sc
               +e->cnt *x[i] *x[k]) /cnt : 0;
      }                         /* (this somewhat complicated formula */
    }                           /* takes missing values into account) */
  }
  if (mode & MVN_CORREL) {      /* correlation coefficients */
    for (i = 0; i < mvn->size; i++) {
      v  = mvn->covs[i];        /* traverse the covariances and */
      c  = mvn->corrs[i] +i;    /* the correlation coefficients */
      *c = 1;                   /* the correlation coefficient of */
      for (k = i; --k >= 0; ) { /* a variable with itself is 1 */
        t    = v[i] *mvn->covs[k][k];
        t    = (t > 0) ? sqrt(t) : 0;
        *--c = (t > 0) ? v[k] /t : 0;
      }                         /* compute the correlation coeffs. */
    }                           /* (set correl. coeff. to 'null' */
  }                             /* if the covariance is null) */
  if ((mode & MVN_DECOM)        /* do Cholesky decomposition */
  &&  (_decom(mvn) != 0)) {     /* and check for success */
    err = -1;                   /* set the error flag */
    mvn->det = pow(EPSILON, mvn->size);
    for (i = mvn->size; --i >= 0; ) {
      c = mvn->decom[i]; c[i] = EPSILON;
      mvn->buf[i] = 1.0/EPSILON;
      for (k = i; ++k < mvn->size; )
        mvn->decom[i][k] = 0;   /* if the decomposition failed, */
    }                           /* set the decomposition of a */
  }                             /* minimal covariance matrix */
  if (mode & MVN_INVERSE)       /* compute the inverse */
    _inv(mvn);                  /* from the Cholesky decomposition */
  return err;                   /* return the error flag */
}  /* mvn_calc() */

/*--------------------------------------------------------------------*/

double mvn_eval (MVNORM *mvn, const double vals[])
{                               /* --- evaluate a normal distribution */
  int    i, k;                  /* loop variables */
  double *d, *x, *p;            /* to traverse vectors and matrices */
  double t;                     /* temporary buffer for exponent */

  assert(mvn && vals);          /* check the function arguments */
  d = mvn->diff +mvn->size;     /* traverse the data values and the */
  x = mvn->exps +mvn->size;     /* corresponding expected values */
  for (vals += (i = mvn->size); --i >= 0; )
    *--d = *--vals - *--x;      /* compute the difference vector d */
  p = mvn->buf +mvn->size;      /* get buffer for intermediate result */
  for (i = mvn->size; --i >= 0; ) {
    *--p = 0;                   /* traverse the matrix columns */
    for (d += k = mvn->size; --k > i; )
      *p += *--d * mvn->inv[k][i];
    for (x = mvn->inv[i] +(k = i+1); --k >= 0; )
      *p += *--d * *--x;        /* calc. product of the diff. vector */
  }                             /* with a column of the inverse */
  d += mvn->size;               /* traverse the diff. vector d again */
  for (t = 0, p += i = mvn->size; --i >= 0; )
    t += *--d * *--p;           /* calc. product with d^T * \Sigma^-1 */
  return mvn->norm *exp(-0.5 *t);      /* return the value */
}  /* mvn_eval() */             /* of the density function */

/*--------------------------------------------------------------------*/

double* mvn_rand (MVNORM *mvn, double drand (void))
{                               /* --- generate random sample point */
  int    i, k;                  /* loop variables */
  double *b, *r;                /* to access the buffer/matrix rows */

  for (b = mvn->buf +(i = mvn->size); --i >= 0; )
    *--b = _normd(drand);       /* generate points from N(0,1)^n */
  for (i = mvn->size; --i >= 0; ) {
    r = mvn->decom[i] +i;       /* traverse the matrix rows */
    b[i] *= *r;                 /* multiply the random vector */
    for (k = i; --k >= 0; )     /* with the Cholesky decomposed */
      b[i] += *--r *b[k];       /* covariance matrix (transposed) */
  }                             /* -> hyperellipsoidal norm. dist. */
  for (b += (i = mvn->size); --i >= 0; )
    *--b += mvn->exps[i];       /* add the expected value vector */
  return b;                     /* return the created sample point */
}  /* mvn_rand() */

/*--------------------------------------------------------------------*/

int mvn_desc (MVNORM *mvn, FILE *file, int offs, int maxlen)
{                               /* --- describe a normal distribution */
  int    i, k;                  /* loop variables */
  double *p;                    /* to traverse the matrix rows */

  assert(mvn && file);          /* check the function arguments */ 

  /* --- expected values --- */
  if (offs > 0) {               /* if the offset is positive, */
    for (k = offs; --k >= 0; )  /* indent the first line */
      fputc(' ', file);         /* to the given offset */
  }                             /* (indent all but the first line */
  offs = abs(offs);             /* if the offset is negative) */
  fputc('[', file);             /* start the expected values vector */
  for (p = mvn->exps, i = 1; i < mvn->size; i++)
    fprintf(file, "%g, ", *p++);/* print the expected values and */
  fprintf(file, "%g],\n", *p);  /* then terminate the vector */

  /* --- (co)variances --- */
  for (i = 0; i < mvn->size; i++) {
    if (i > 0)                  /* if this is not the first row, */
      fputs(",\n", file);       /* print a comma and start new line */
    for (k = offs; --k >= 0; )  /* indent the line */
      fputc(' ', file);         /* to the given offset */
    fputc('[', file);           /* start a new matrix row */
    for (p = mvn->covs[i], k = 0; k < i; k++)
      fprintf(file, "%g, ", *p++);      /* print the covariances */
    fprintf(file, "%g]", *p);   /* and then terminate the vector */
  }                             /* by printing the variance */

  return ferror(file) ? -1 : 0; /* return file write status */
}  /* mvn_desc() */

/*--------------------------------------------------------------------*/
#ifdef MVN_PARSE

static int _get_exps (MVNORM *mvn, SCAN *scan)
{                               /* --- parse expected values */
  int    i;                     /* loop variable */
  MVNROW *row;                  /* to traverse the matrix rows */
  double t;                     /* temporary buffer */

  assert(mvn && scan);          /* check the function arguments */
  GET_CHR('[');                 /* consume '(' */
  for (i = 0; i < mvn->size; i++) {
    if (i > 0) {GET_CHR(',');}  /* check for a comma and consume it */
    if (sc_token(scan) != T_NUM) ERROR(E_NUMEXP);
    mvn->exps[i] = t = atof(sc_value(scan));
    GET_TOK();                  /* get and consume expected value */
    row = mvn->rows[i];         /* recompute the cum of values */
    row->sv = t *row->cnt;      /* from the expected value */
  }                             /* and the number of cases */
  GET_CHR(']');                 /* consume ')' */
  GET_CHR(',');                 /* consume ',' */
  return 0;                     /* return 'ok' */
}  /* _get_exps() */

/*--------------------------------------------------------------------*/

static int _get_covs (MVNORM *mvn, SCAN *scan)
{                               /* --- parse (co)variances */
  int     i, k;                 /* loop variables */
  MVNROW  *row;                 /* to traverse the matrix rows */
  MVNELEM *e;                   /* to traverse the matrix elements */
  double  *v;                   /* to traverse the (co)variances */

  assert(mvn && scan);          /* check the function arguments */
  for (i = 0; i < mvn->size; i++) {
    row = mvn->rows[i];         /* traverse the matrix rows */
    if (i > 0) {GET_CHR(',');}  /* check for a comma between rows */
    GET_CHR('[');               /* consume '(' */
    e = row->elems;             /* traverse the statistics and */
    v = mvn->covs[i];           /* the covariances of each row */
    for (k = 0; k < i; e++, v++, k++) {
      if (sc_token(scan) != T_NUM) ERROR(E_NUMEXP);
      *v = atof(sc_value(scan));/* get the covariance */
      e->src = *v *(row->cnt -1) +row->sv *mvn->exps[k];
      GET_TOK();                /* recompute mixed sum, */
      GET_CHR(',');             /* consume covariance, */
    }                           /* and consume ',' */
    if (sc_token(scan) != T_NUM) ERROR(E_NUMEXP);
    *v = atof(sc_value(scan));  /* get and check the variance */
    if (*v < 0)                  ERROR(E_NUMBER);
    GET_TOK();                  /* consume variance */
    GET_CHR(']');               /* consume ')' */
    row->sv2 = *v *(row->cnt -1) +row->sv *mvn->exps[i];
    if (row->sv2 < 0) row->sv2 = 0;
    for (k = i; --k >= 0; ) {   /* recompute sum of squares and */
      --e;                      /* traverse the elements again */
      e->cnt = row->cnt;        /* set the number of cases */
      e->sr  = row->sv;  e->sc  = mvn->rows[k]->sv;
      e->sr2 = row->sv2; e->sc2 = mvn->rows[k]->sv2;
    }                           /* get the sums of values and */
  }                             /* the sums of squared values */
  return 0;                     /* return 'ok' */
}  /* _get_covs() */

/*--------------------------------------------------------------------*/

static int _get_poss (MVNORM *mvn, SCAN *scan)
{                               /* --- get possibilistic parameter */
  assert(mvn && scan);          /* check the function arguments */
  if (sc_token(scan) != ',') {  /* if no comma follows, */
    mvn->poss = 1; return 0; }  /* set a default value */
  GET_TOK();                    /* consume ',' */
  if (sc_token(scan) != T_NUM) ERROR(E_NUMEXP);
  mvn->poss = atof(sc_value(scan));
  if (mvn->poss <= 0)          ERROR(E_NUMBER);
  GET_TOK();                    /* get and consume the parameter */
  return 0;                     /* return 'ok' */
}  /* _get_poss() */

/*--------------------------------------------------------------------*/

int mvn_parse (MVNORM *mvn, SCAN *scan, double cnt)
{                               /* --- parse a normal distribution */
  int    i;                     /* loop variable, function result */
  MVNROW **row;                 /* to traverse the matrix rows */

  assert(mvn && scan);          /* check the function arguments */
  if (cnt < 2) cnt = 2;         /* adapt the number of cases */
  for (row = mvn->rows +(i = mvn->size); --i >= 0; )
    (*--row)->cnt = cnt;        /* set number of cases */
  pa_init(scan);                /* initialize parsing */
  i = _get_exps(mvn, scan);     /* read the expected values */
  if (i < 0) return i;
  i = _get_covs(mvn, scan);     /* read the (co)variances */
  if (i < 0) return i;
  i = _get_poss(mvn, scan);     /* read the possibilistic parameter */
  if (i < 0) return i;
  return 0;                     /* return 'ok' */
}  /* mvn_parse() */

#endif