randvar.c

来自「General Hidden Markov Model Library 一个通用」· C语言 代码 · 共 942 行 · 第 1/2 页

/*******************************************************************************
*
*       This file is part of the General Hidden Markov Model Library,
*       GHMM version 0.8_beta1, see http://ghmm.org
*
*       Filename: ghmm/ghmm/randvar.c
*       Authors:  Bernhard Knab, Benjamin Rich, Janne Grunau
*
*       Copyright (C) 1998-2004 Alexander Schliep
*       Copyright (C) 1998-2001 ZAIK/ZPR, Universitaet zu Koeln
*       Copyright (C) 2002-2004 Max-Planck-Institut fuer Molekulare Genetik,
*                               Berlin
*
*       Contact: [e-mail]
*
*       This library is free software; you can redistribute it and/or
*       modify it under the terms of the GNU Library General Public
*       License as published by the Free Software Foundation; either
*       version 2 of the License, or (at your option) any later version.
*
*       This library is distributed in the hope that it will be useful,
*       but WITHOUT ANY WARRANTY; without even the implied warranty of
*       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
*       Library General Public License for more details.
*
*       You should have received a copy of the GNU Library General Public
*       License along with this library; if not, write to the Free
*       Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*
*       This file is version $Revision: 1650 $
*                       from $Date: 2006-07-20 17:57:27 +0200 (Thu, 20 Jul 2006) $
*             last change by $Author: grunau $.
*
*******************************************************************************/

#ifdef WIN32
#  include "win_config.h"
#endif

#ifdef HAVE_CONFIG_H
#  include "../config.h"
#endif

#include <math.h>
#include <float.h>
#include <string.h>
#ifdef HAVE_LIBPTHREAD
# include <pthread.h>
#endif /* HAVE_LIBPTHREAD */

#include "ghmm.h"
#include "mes.h"
#include "mprintf.h"
#include "randvar.h"
#include "rng.h"

#ifdef DO_WITH_GSL

# include <gsl/gsl_math.h>
# include <gsl/gsl_sf_erf.h>
# include <gsl/gsl_randist.h>

#endif /* DO_WITH_GSL */


#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)
#else
static double ighmm_erf (double x);
static double ighmm_erfc (double x);
#endif /* check for ISO C99 */


/* A list of already calculated values of the density function of a 
   N(0,1)-distribution, with x in [0.00, 19.99] */
#define PDFLEN 2000
#define X_STEP_PDF 0.01         /* step size */
#define X_FAKT_PDF 100          /* equivalent to step size */
static double pdf_stdnormal[PDFLEN];
static int pdf_stdnormal_exists = 0;

/* A list of already calulated values PHI of the Gauss distribution is
   read in, x in [-9.999, 0] */
#define X_STEP_PHI 0.001        /* step size */
#define X_FAKT_PHI 1000         /* equivalent to step size */
static int PHI_len = 0;
static double x_PHI_1 = -1.0;

#ifndef M_SQRT1_2
#define M_SQRT1_2  0.70710678118654752440084436210
#endif


/*============================================================================*/
/* needed by ighmm_gtail_pmue_interpol */

double ighmm_rand_get_xfaktphi ()
{
  return X_FAKT_PHI;
}

double ighmm_rand_get_xstepphi ()
{
  return X_STEP_PHI;
}

double ighmm_rand_get_philen ()
{
#ifdef DO_WITH_GSL
  return PHI_len;
#else
  return ighmm_rand_get_xPHIless1 () / X_STEP_PHI;
#endif
}


/*============================================================================*/
double ighmm_rand_get_PHI (double x)
{
#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)
  return (erf (x * M_SQRT1_2) + 1.0) / 2.0;
#else
  return (ighmm_erf (x * M_SQRT1_2) + 1.0) / 2.0;
#endif
}                               /* randvar_get_PHI */


/*============================================================================*/
/* When is PHI[x,0,1] == 1? */
double ighmm_rand_get_xPHIless1 ()
{
# define CUR_PROC "ighmm_rand_get_xPHIless1"

  if (x_PHI_1 == -1) {
    double low, up, half;
    low = 0;
    up = 100;
    while (up - low > 0.001) {
      half = (low + up) / 2.0;
      if (ighmm_rand_get_PHI (half) < 1.0)
        low = half;
      else
        up = half;
    }
    x_PHI_1 = low;
  }
  return (x_PHI_1);

# undef CUR_PROC
}

/*============================================================================*/
double ighmm_rand_get_1overa (double x, double mean, double u)
{
  /* Calulates 1/a(x, mean, u), with a = the integral from x til \infty over
     the Gauss density function */
# define CUR_PROC "ighmm_rand_get_1overa"

  double erfc_value;

  if (u <= 0.0) {
    GHMM_LOG(LCONVERTED, "u <= 0.0 not allowed\n");
    goto STOP;
  }

#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)
  erfc_value = erfc ((x - mean) / sqrt (u * 2));
#else
  erfc_value = ighmm_erfc ((x - mean) / sqrt (u * 2));
#endif

  if (erfc_value <= DBL_MIN) {
    ighmm_mes (MES_WIN, "a ~= 0.0 critical! (mue = %.2f, u =%.2f)\n", mean, u);
    return (erfc_value);
  }
  else
    return (2.0 / erfc_value);

STOP:
  return (-1.0);
# undef CUR_PROC
}                               /* ighmm_rand_get_1overa */


/*============================================================================*/
/* REMARK:
   The calulation of this density function was testet, by calculating the 
   following integral sum for arbitrary mue and u:
     for (x = 0, x < ..., x += step(=0.01/0.001/0.0001)) 
       isum += step * ighmm_rand_normal_density_pos(x, mue, u);
   In each case, the sum "converged" evidently towards 1!
   (BK, 14.6.99)
   CHANGE:
   Truncate at -EPS_NDT (const.h), so that x = 0 doesn't lead to a problem.
   (BK, 15.3.2000)
*/
double ighmm_rand_normal_density_pos (double x, double mean, double u)
{
# define CUR_PROC "ighmm_rand_normal_density_pos"
  return ighmm_rand_normal_density_trunc (x, mean, u, -GHMM_EPS_NDT);
# undef CUR_PROC
}                               /* double ighmm_rand_normal_density_pos */


/*============================================================================*/
double ighmm_rand_normal_density_trunc (double x, double mean, double u,
                                     double a)
{
# define CUR_PROC "ighmm_rand_normal_density_trunc"
#ifndef DO_WITH_GSL
  double c;
#endif /* DO_WITH_GSL */

  if (u <= 0.0) {
    GHMM_LOG(LCONVERTED, "u <= 0.0 not allowed\n");
    goto STOP;
  }
  if (x < a)
    return (0.0);

#ifdef DO_WITH_GSL
  /* move mean to the right position */
  /* double gsl_ran_gaussian_tail_pdf (double x, double a, double sigma) */
  return gsl_ran_gaussian_tail_pdf (x - mean, a - mean, sqrt (u));
#else
  if ((c = ighmm_rand_get_1overa (a, mean, u)) == -1) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  };
  return (c * ighmm_rand_normal_density (x, mean, u));
#endif /* DO_WITH_GSL */

STOP:
  return (-1.0);
# undef CUR_PROC
}                               /* double ighmm_rand_normal_density_trunc */


/*============================================================================*/
double ighmm_rand_normal_density (double x, double mean, double u)
{
# define CUR_PROC "ighmm_rand_normal_density"
#ifndef DO_WITH_GSL
  double expo;
#endif
  if (u <= 0.0) {
    GHMM_LOG(LCONVERTED, "u <= 0.0 not allowed\n");
    goto STOP;
  }
  /* The denominator is possibly < EPS??? Check that ? */
#ifdef DO_WITH_GSL
  /* double gsl_ran_gaussian_pdf (double x, double sigma) */
  return gsl_ran_gaussian_pdf (x - mean, sqrt (u));
#else
  expo = exp (-1 * m_sqr (mean - x) / (2 * u));
  return (1 / (sqrt (2 * PI * u)) * expo);
#endif

STOP:
  return (-1.0);
# undef CUR_PROC
}                               /* double ighmm_rand_normal_density */


/*============================================================================*/
double ighmm_rand_uniform_density (double x, double max, double min)
{
# define CUR_PROC "ighmm_rand_uniform_density"
  double prob;
  if (max <= min) {
    GHMM_LOG(LCONVERTED, "max <= min not allowed \n");
    goto STOP;
  }
  prob = 1.0/(max-min);

  if ( (x <= max) && (x>=min) ){
    return prob;
  }else{
    return 0.0;
  }
STOP:
  return (-1.0);
# undef CUR_PROC
}                               /* double ighmm_rand_uniform_density */


/*============================================================================*/
/* special ghmm_cmodel pdf need it: smo->density==normal_approx: */
/* generates a table of of aequidistant samples of gaussian pdf */

static int randvar_init_pdf_stdnormal ()
{
# define CUR_PROC "randvar_init_pdf_stdnormal"
  int i;
  double x = 0.00;
  for (i = 0; i < PDFLEN; i++) {
    pdf_stdnormal[i] = 1 / (sqrt (2 * PI)) * exp (-1 * x * x / 2);
    x += (double) X_STEP_PDF;
  }
  pdf_stdnormal_exists = 1;
  /* printf("pdf_stdnormal_exists = %d\n", pdf_stdnormal_exists); */
  return (0);
# undef CUR_PROC
}                               /* randvar_init_pdf_stdnormal */


double ighmm_rand_normal_density_approx (double x, double mean, double u)
{
# define CUR_PROC "ighmm_rand_normal_density_approx"
#ifdef HAVE_LIBPTHREAD
  static pthread_mutex_t lock;
#endif /* HAVE_LIBPTHREAD */
  int i;
  double y, z, pdf_x;
  if (u <= 0.0) {
    GHMM_LOG(LCONVERTED, "u <= 0.0 not allowed\n");
    goto STOP;
  }
  if (!pdf_stdnormal_exists) {
#ifdef HAVE_LIBPTHREAD
    pthread_mutex_lock (&lock); /* Put on a lock, because the clustering is parallel   */
#endif /* HAVE_LIBPTHREAD */
    randvar_init_pdf_stdnormal ();
#ifdef HAVE_LIBPTHREAD
    pthread_mutex_unlock (&lock);       /* Take the lock off */
#endif /* HAVE_LIBPTHREAD */
  }
  y = 1 / sqrt (u);
  z = fabs ((x - mean) * y);
  i = (int) (z * X_FAKT_PDF);
  /* linear interpolation: */
  if (i >= PDFLEN - 1) {
    i = PDFLEN - 1;
    pdf_x = y * pdf_stdnormal[i];
  }
  else
    pdf_x = y * (pdf_stdnormal[i] +
                 (z - i * X_STEP_PDF) *
                 (pdf_stdnormal[i + 1] - pdf_stdnormal[i]) / X_STEP_PDF);
  return (pdf_x);
STOP:
  return (-1.0);
# undef CUR_PROC
}                               /* double ighmm_rand_normal_density_approx */


/*============================================================================*/
double ighmm_rand_std_normal (int seed)
{
# define CUR_PROC "ighmm_rand_std_normal"
  if (seed != 0) {
    GHMM_RNG_SET (RNG, seed);
    return (1.0);
  }
  else {
#ifdef DO_WITH_GSL
    return (gsl_ran_gaussian (RNG, 1.0));
#else
    /* Use the polar Box-Mueller transform */
    /*
       double x, y, r2;

       do {
       x = 2.0 * GHMM_RNG_UNIFORM(RNG) - 1.0;
       y = 2.0 * GHMM_RNG_UNIFORM(RNG) - 1.0;
       r2 = (x * x) + (y * y);
       } while (r2 >= 1.0);

       return x * sqrt((-2.0 * log(r2)) / r2);
     */

    double r2, theta;

    r2 = -2.0 * log (GHMM_RNG_UNIFORM (RNG));   /* r2 ~ chi-square(2) */
    theta = 2.0 * PI * GHMM_RNG_UNIFORM (RNG);  /* theta ~ uniform(0, 2 \pi) */
    return sqrt (r2) * cos (theta);
#endif
  }
# undef CUR_PROC
}                               /* ighmm_rand_std_normal */


/*============================================================================*/
double ighmm_rand_normal (double mue, double u, int seed)
{
# define CUR_PROC "ighmm_rand_normal"
  if (seed != 0) {
    GHMM_RNG_SET (RNG, seed);
    return (1.0 * sqrt (u) + mue);
  }
  else {
#ifdef DO_WITH_GSL
    return (gsl_ran_gaussian (RNG, sqrt (u)) + mue);
#else
    double x;
    x = sqrt (u) * ighmm_rand_std_normal (seed) + mue;
    return (x);
#endif
  }
# undef CUR_PROC
}                               /* ighmm_rand_normal */


/*============================================================================*/
#ifndef DO_WITH_GSL
# define C0 2.515517
# define C1 0.802853
# define C2 0.010328
# define D1 1.432788
# define D2 0.189269
# define D3 0.001308
#endif

double ighmm_rand_normal_right (double a, double mue, double u, int seed)
{
# define CUR_PROC "ighmm_rand_normal_right"
  double x = -1;
  double sigma;
#ifdef DO_WITH_GSL
  double s;
#else
  double U, Us, Us1, Feps, Feps1, t, T;
#endif

  if (u <= 0.0) {
    GHMM_LOG(LCONVERTED, "u <= 0.0 not allowed\n");
    goto STOP;
  }
  sigma = sqrt (u);

  if (seed != 0) {
    GHMM_RNG_SET (RNG, seed);
    return (1.0);
  }

#ifdef DO_WITH_GSL
  /* up to version 0.8 gsl_ran_gaussian_tail can not handle negative cutoff */
#define GSL_RAN_GAUSSIAN_TAIL_BUG 1
#ifdef GSL_RAN_GAUSSIAN_TAIL_BUG
  s = (-mue) / sigma;
  if (s < 1) {
    do {
      x = gsl_ran_gaussian (RNG, 1.0);
    }
    while (x < s);
    return x * sigma + mue;
  }
#endif /* GSL_RAN_GAUSSIAN_TAIL_BUG */
  /* move boundary to lower values in order to achieve maximum at mue
     gsl_ran_gaussian_tail(generator, lower_boundary, sigma)
   */
  return gsl_ran_gaussian_tail (RNG, a - mue, sqrt (u)) + mue;

#else /* DO_WITH_GSL */
  /* Method: Generate Gauss-distributed random nunbers (with GSL-lib.),
     until a positive one is found -> not very effective if mue << 0
     while (x < 0.0) {
     x = sigma * ighmm_rand_std_normal(seed) + mue;
     } */

  /* Inverse transformation with restricted sampling by Fishman */
  U = GHMM_RNG_UNIFORM (RNG);
  Feps = ighmm_rand_get_PHI (-(a + mue) / sigma);
  Us = Feps + (1 - Feps) * U;
  /* Numerically better: 1-Us = 1-Feps - (1-Feps)*U, therefore: 
     Feps1 = 1-Feps, Us1 = 1-Us */
  Feps1 = ighmm_rand_get_PHI ((a + mue) / sigma);
  Us1 = Feps1 - Feps1 * U;
  t = m_min (Us, Us1);
  t = sqrt (-log (t * t));
  T =