mvnorm.c

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

/*----------------------------------------------------------------------
  File    : mvnorm.c
  Contents: Multivariate normal distribution estimation and management
  Author  : Christian Borgelt
  History : 2000.11.10 file created from files correl.c and cluster.c
            2000.11.26 first version completed
            2001.05.24 possibilistic parameter added
            2001.05.26 handling of roundoff errors improved
            2001.07.15 some assertions added
            2004.04.15 treatment of decomposition failure improved
            2004.08.12 adapted to new module parse
            2005.09.05 bug in function _decom (recomputation) fixed
----------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#ifdef MVN_PARSE
#include <string.h>
#endif
#include <math.h>
#include <assert.h>
#include "mvnorm.h"
#ifdef STORAGE
#include "storage.h"
#endif

/*----------------------------------------------------------------------
  Preprocessor Definitions
----------------------------------------------------------------------*/
#define	M_PI        3.14159265358979323846  /* \pi */
#define EPSILON     1e-12       /* to handle roundoff errors */

/*----------------------------------------------------------------------
  Auxiliary Functions
----------------------------------------------------------------------*/

static double _normd (double drand (void))
{                               /* --- compute N(0,1) distrib. number */
  static double b;              /* buffer for random number */
  double x, y, r;               /* coordinates and radius */

  if (b != 0.0) {               /* if the buffer is full, */
    x = b; b = 0; return x; }   /* return the buffered number */
  do {                          /* pick a random point */
    x = 2.0*drand()-1.0;        /* in the unit square [-1,1]^2 */
    y = 2.0*drand()-1.0;        /* and check whether it lies */
    r = x*x +y*y;               /* inside the unit circle */
  } while ((r > 1) || (r == 0));
  r = sqrt(-2*log(r)/r);        /* factor for Box-Muller transform */
  b = x *r;                     /* save one of the random numbers */
  return y *r;                  /* and return the other */
}  /* _normd() */

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

static int _decom (MVNORM *mvn)
{                               /* --- Cholesky decomposition */
  int    i, j, k;               /* loop variables */
  double *s, *r;                /* to traverse the matrix rows */
  double t, q;                  /* temporary buffers */

  assert(mvn);                  /* check the function argument */
  for (i = 0; i < mvn->size; i++) {
    s = mvn->decom[i];          /* traverse the matrix rows */
    for (t = mvn->covs[i][k = i]; --k >= 0; )
      t -= s[k] *s[k];          /* compute square of diag. element */
    if (t <= 0) return -1;      /* matrix is not positive definite */
    s[i] = t = sqrt(t);         /* store the diagonal element */
    if (t <= 0) return -1;      /* matrix is not positive definite */
    mvn->buf[i] = q = 1/t;      /* compute the corresponding factor */
    for (j = i; ++j < mvn->size; ) {
      r = mvn->decom[j];        /* traverse the row's other columns */
      for (t = mvn->covs[j][k = i]; --k >= 0; )
        t -= s[k] *r[k];        /* compute the off-diagonal elements */
      r[i] = t *q;              /* column by column */
    }                           /* (-> lower triangular matrix) */
  }
  for (t = mvn->decom[0][0], i = mvn->size; --i > 0; )
    t *= mvn->decom[i][i];      /* multiply the diagonal elements */
  mvn->det = t*t;               /* to compute the determinant */
  if (mvn->det <= 0) return -1; /* check for a positive determinant */
  return 0;                     /* return 'ok' */
}  /* _decom() */

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

static void _inv (MVNORM *mvn)
{                               /* --- invert matrix from Cholesky d. */
  int    i, k, n;               /* loop variables */
  double *d;                    /* reciprocals of diagonal elements */
  double *s, *r;                /* to traverse the matrix elements */
  double t;                     /* temporary buffer */

  assert(mvn);                  /* check the function argument */
  d = mvn->buf;                 /* get reciprocals of diag. elements */
  for (n = mvn->size; --n >= 0; ) {
    r = mvn->diff;              /* traverse the unit vectors */
    r[n] = d[n];                /* and compute their originals */
    for (i = n; ++i < mvn->size; ) {
      s = mvn->decom[i];        /* traverse the matrix rows */
      for (t = 0, k = i; --k >= n; )
        t -= r[k] *s[k];        /* do forward substitution with */
      r[i] = t *d[i];           /* the lower triangular matrix L */
    }                           /* (store result in r = mvn->diff) */
    for (i = mvn->size; --i >= n; ) {
      for (t = r[i], k = mvn->size; --k > i; )
        t -= mvn->decom[k][i] *mvn->inv[k][n];
      mvn->inv[i][n] = t*d[i];  /* do backward substitution with */
    }                           /* the upper triangular matrix L^T */
  }                             /* (using only the lower triangle) */
  mvn->norm = 1/sqrt(mvn->det *pow(2*M_PI, mvn->size));
}  /* _inv() */                 /* compute normalization factor */

/*----------------------------------------------------------------------
  Main Functions
----------------------------------------------------------------------*/

static MVNORM* _create (int size)
{                               /* --- create a normal distribution */
  int    i;                     /* loop variable */
  MVNORM *mvn;                  /* created normal distribution */
  MVNROW **row;                 /* to traverse the matrix rows */
  double *p;                    /* to traverse the matrix rows */

  if (size <= 0) size = 1;      /* check and adapt the matrix size */
  mvn = (MVNORM*)calloc(1, sizeof(MVNORM) +(size-1) *sizeof(MVNROW*));
  if (!mvn) return NULL;        /* create the base structure */
  mvn->size  = size;            /* and store the matrix size */
  mvn->covs  = (double**)malloc(4*size *sizeof(double*));
  if (!mvn->covs) { free(mvn); return NULL; }
  mvn->corrs = mvn->covs  +size;   /* create vectors of pointers */
  mvn->decom = mvn->corrs +size;   /* for matrix rows */
  mvn->inv   = mvn->decom +size;
  mvn->exps  = (double*)malloc((2*size*size +5*size) *sizeof(double));
  if (!mvn->exps) { mvn_delete(mvn); return NULL; }
  p = mvn->exps +size;          /* create and organize data vectors */
  for (i = size; --i >= 0; ) { mvn->covs[i]  = p; p += i+1; }
  for (i = size; --i >= 0; ) { mvn->corrs[i] = p; p += i+1; }
  for (i = size; --i >= 0; ) { mvn->decom[i] = p; p += i+1; }
  for (i = size; --i >= 0; ) { mvn->inv[i]   = p; p += i+1; }
  mvn->diff = p;                /* set buffers for difference from */
  mvn->buf  = p +size;          /* center and intermediate results */
  for (row = mvn->rows +size; --size >= 0; ) {
    i = (size > 0) ? size-1 : 0;
    *--row = (MVNROW*)malloc(sizeof(MVNROW) +i *sizeof(MVNELEM));
    if (!*row) { mvn_delete(mvn); return NULL; }
  }                             /* create the statistics rows */
  return mvn;                   /* and return it */
}  /* _create() */

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

MVNORM* mvn_create (int size)
{                                /* --- create a normal distribution */
  MVNORM *mvn;                   /* created normal distribution */

  mvn = _create(size);           /* create a normal distribution */
  if (!mvn) return NULL;
  mvn_clear(mvn);                /* clear the created normal dist. */
  return mvn;                    /* and return it */
}  /* mvn_create() */

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

MVNORM* mvn_clone (const MVNORM *mvn)
{                               /* --- clone a normal distribution */
  int     i, k;                 /* loop variables */
  MVNORM  *clone;               /* created clone */
  MVNROW  *dr; const MVNROW  *sr;  /* to traverse the matrix rows */
  MVNELEM *de; const MVNELEM *se;  /* to traverse the matrix elements */
  double  *d;  const double  *s;   /* to traverse the parameters */

  assert(mvn);                  /* check the function argument */
  clone = _create(mvn->size);   /* create a normal distribution */
  if (!clone) return NULL;      /* of the same size */
  clone->det  = mvn->det;       /* copy the determinant */
  clone->norm = mvn->norm;      /* and the normalization factor */
  for (i = mvn->size; --i >= 0; ) {
    sr = mvn->rows[i];          /* traverse the matrix rows */
    dr = clone->rows[i];        /* of source and destination */
    *dr = *sr;                  /* and copy them */
    se = sr->elems +i;          /* traverse the row elements */
    de = dr->elems +i;          /* and copy them */
    for (k = i; --k >= 0; ) *--de = *--se;
  }
  i = 2*mvn->size*mvn->size +5*mvn->size;
  s = mvn->exps +i;             /* traverse and copy the vectors of */
  d = clone->exps +i;           /* estimated parameters, the Cholesky */
  while (--i >= 0) *--d = *--s; /* decomposition, and the inverse */
  return clone;                 /* return the created clone */
}  /* mvn_clone() */

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

void mvn_delete (MVNORM *mvn)
{                               /* --- delete a normal distribution */
  int    i;                     /* loop variable */
  MVNROW **row;                 /* to traverse the matrix rows */

  assert(mvn);                  /* check the function argument */
  for (row = mvn->rows +(i = mvn->size); --i >= 0; )
    if (*--row) free(*row);     /* delete all statistics rows, */
  if (mvn->exps) free(mvn->exps);    /* the parameter vectors, */
  if (mvn->covs) free(mvn->covs);    /* the pointer vectors, */
  free(mvn);                         /* and the base structure */
}  /* mvn_delete() */

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

void mvn_clear (MVNORM *mvn)
{                               /* --- clear a normal distribution */
  int     i, k;                 /* loop variables */
  MVNROW  *row;                 /* to traverse the matrix rows */
  MVNELEM *e;                   /* to traverse the matrix elements */

  assert(mvn);                  /* check the function argument */
  for (i = mvn->size; --i >= 0; ) {
    row = mvn->rows[i];         /* traverse the matrix rows and */
    row->cnt = row->sv = row->sv2 = 0;        /* clear the sums */
    for (e = row->elems +(k = i); --k >= 0; ) {
      (--e)->cnt = 0; e->sr = e->sr2 = e->sc = e->sc2 = e->src = 0; }
  }                             /* clear all sums of the */
}  /* mvn_clear() */            /* matrix elements */

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

void mvn_add (MVNORM *mvn, const double vals[], double cnt)
{                               /* --- add a value vector */
  int     i, k;                 /* loop variables */
  MVNROW  *r1, *r2;             /* to traverse the matrix rows */
  MVNELEM *e;                   /* to traverse the matrix elements */

  assert(mvn && vals && (cnt >= 0)); /* check the function arguments */
  for (i = 0; i < mvn->size; i++) {
    if (vals[i] <= MVN_NULL)    /* traverse the matrix rows, */
      continue;                 /* but skip those rows, for */
    r1 = mvn->rows[i];          /* which the value is null */
    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] <= MVN_NULL)  /* traverse the preceding rows */
        continue;               /* but skip those rows, for */
      r2 = mvn->rows[k];        /* which the value is null */
      (--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_add() */

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

void mvn_addx (MVNORM *mvn, const double vals[], double cnt)
{                               /* --- add a value vector */
  int     i, k;                 /* loop variables */
  MVNROW  *r1, *r2;             /* to traverse the matrix rows */
  MVNELEM *e;                   /* to traverse the matrix elements */

  assert(mvn && vals && (cnt >= 0)); /* check the function arguments */
  for (i = 0; i < mvn->size; i++) {
    if (vals[i] <= 0)           /* traverse the matrix rows, */
      continue;                 /* but skip those rows, for */
    r1 = mvn->rows[i];          /* which the value is not positive */