bvnorm.c

BCView - Bayes Classifier Visualization Download xbcview Linux executable (218 kb) wbcview.exe W

/*----------------------------------------------------------------------
  File    : bvnorm.c
  Contents: Bivariate normal distribution management
  Author  : Christian Borgelt
  History : 14.01.2000 file created
            31.01.2000 function bvn_ellphi added
            01.02.2000 functions bvn_ellx2y and bvn_elly2x added
            06.03.2001 less restrictive parameter checking in bvn_init
            27.04.2004 function bvn_dist added
----------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <assert.h>
#include "bvnorm.h"

/*----------------------------------------------------------------------
  Preprocessor Definitions
----------------------------------------------------------------------*/
#define	M_PI  3.14159265358979323846  /* \pi */

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

int bvn_init (BVNORM *bvn, double prob, double ex, double ey,
              double cxx, double cyy, double cxy)
{                               /* --- create a normal distribution */
  double t;                     /* temporary buffer */

  assert(bvn && (prob >= 0));   /* check the function arguments */
  if ((cxx < 0) || (cyy < 0))   /* if either variance is negative, */
    return -1;                  /* abort the function */
  t = cxx*cyy -cxy*cxy;         /* compute and check the determinant */
  if (t < 0) return -1;         /* of the covariance matrix */
  bvn->det  =  t;               /* and store it */
  bvn->prob =  prob;            /* note prior probability */
  bvn->ex   =  ex;              /* and expected values */
  bvn->ey   =  ey;
  bvn->dx   =  sqrt(cxx);       /* compute and store */
  bvn->dy   =  sqrt(cyy);       /* standard deviations */
  bvn->cxx  =  cxx;             /* store the elements of */
  bvn->cyy  =  cyy;             /* the covariance matrix */
  bvn->cxy  =  cxy;             
  if (t <= 0)                   /* if the determinant is zero, */
    bvn->ixx = bvn->iyy = bvn->ixy = 0;   /* clear the inverse */
  else {                        /* otherwise */
    bvn->ixx =  cyy/t;          /* compute and store the */
    bvn->iyy =  cxx/t;          /* elements of the inverse */
    bvn->ixy = -cxy/t;          /* of the covariance matrix */
  }
  t = sqrt(t);                  /* compute the normalization factor */
  bvn->norm =  prob/((t > 0) ? 2 *M_PI *t : DBL_MIN);
  if (bvn->dx <= 0)             /* if the x-variance is zero, clear */
    bvn->dxx = bvn->dxy = bvn->dyy = 0;        /* the decomposition */
  else {                        /* otherwise */
    bvn->dxx = bvn->dx;         /* compute a Cholesky decomposition */
    bvn->dxy = cxy/bvn->dxx;    /* of the covariance matrix */
    bvn->dyy = t  /bvn->dxx;    /* (for bvn_ellx2y and bvn_elly2x) */
  }
  return 0;                     /* return `ok' */
}  /* bvn_init() */

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

double bvn_dist (BVNORM *bvn, double x, double y)
{                               /* --- squared Mahalanobis distance */
  register double t;            /* temporary buffer */

  assert(bvn);                  /* check the function argument */
  x -= bvn->ex; y -= bvn->ey;   /* shift to expected value */
  t = bvn->ixx *(x*x) +bvn->iyy *(y*y);
  if (bvn->ixy != 0) t += 2 *bvn->ixy *(x*y);
  return t;                     /* compute the Mahalanobis distance */
}  /* bvn_dist() */

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

double bvn_eval (BVNORM *bvn, double x, double y)
{                               /* --- evaluate a normal distribution */
  register double t;            /* temporary buffer */

  assert(bvn);                  /* check the function argument */
  x -= bvn->ex; y -= bvn->ey;   /* shift to expected value */
  t = bvn->ixx *(x*x) +bvn->iyy *(y*y);
  if (bvn->ixy != 0) t += 2 *bvn->ixy *(x*y);
  return bvn->norm *exp(-0.5*t);/* compute the function value */
}  /* bvn_eval() */

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

void bvn_ellx2y (BVNORM *bvn, double k, double x, double *y1,double *y2)
{                               /* --- compute points on k sigma ell. */
  double t, u, d;               /* temporary buffers */

  assert(bvn && y1 && y2);      /* check the function arguments */
  if (bvn->cxx <= 0) { *y1 = *y2 = bvn->ey; return; }
  u = bvn->cxy *(d = x-bvn->ex);/* compute base value */
  t = bvn->det *(k*k *bvn->cxx -d*d);
  t = (t > 0) ? sqrt(t) : 0;    /* compute offset (+/- part) */
  *y1 = (u-t)/bvn->cxx +bvn->ey;/* compute y-coordinates */
  *y2 = (u+t)/bvn->cxx +bvn->ey;/* of ellipse points at x */
}  /* bvn_ellx2y() */

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

void bvn_elly2x (BVNORM *bvn, double k, double y, double *x1,double *x2)
{                               /* --- compute points on k sigma ell. */
  double t, u, d;               /* temporary buffers */

  assert(bvn && x1 && x2);      /* check the function arguments */
  if (bvn->cyy <= 0) { *x1 = *x2 = bvn->ex; return; }
  u = bvn->cxy *(d = y-bvn->ey);/* compute base value */
  t = bvn->det *(k*k *bvn->cyy -d*d);
  t = (t > 0) ? sqrt(t) : 0;    /* compute offset (+/- part) */
  *x1 = (u-t)/bvn->cyy +bvn->ex;/* compute x-coordinates */
  *x2 = (u+t)/bvn->cyy +bvn->ex;/* of ellipse points at y */
}  /* bvn_elly2x() */

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

void bvn_ellphi (BVNORM *bvn, double k, double phi, double *x,double *y)
{                               /* --- compute point on k sigma ell. */
  double tx, ty;                /* temporary buffers */

  assert(bvn && x && y);        /* check the function arguments */
  tx = k *cos(phi);             /* compute coords. of point on circle */
  ty = k *sin(phi);             /* transform it with the inv. matrix */
  *x = tx *bvn->dxx               +bvn->ex;
  *y = ty *bvn->dxy +ty *bvn->dyy +bvn->ey;
}  /* bvn_ellphi() */