randvar.c

一个通用的隐性马尔可夫C代码库 开发环境:C语言 简要说明:这是一个通用的隐性马尔可夫C代码库

    a = 1 - phi_z;
    if (a > DBL_MIN) a = 1 / a;
    else {
      a = 0.0;
      mes(MES_WIN, "a ~= 0.0 critical! (mue = %.2f, u =%.2f)\n",
	  mean, u); /* goto STOP; */
    }
  }
  return a;
#endif

STOP:
  return(-1.0);
# undef CUR_PROC
} /* randvar_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 * randvar_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 randvar_normal_density_pos(double x, double mean, double u) { 
# define CUR_PROC "randvar_normal_density_pos"
  return randvar_normal_density_trunc(x, mean, u, -EPS_NDT);
# undef CUR_PROC
} /* double randvar_normal_density_pos */


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

  if (u <= 0.0) {
    mes_prot("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 = randvar_get_1overa(a, mean, u)) == -1) 
    {mes_proc(); goto STOP;};  
  return(c * randvar_normal_density(x, mean, u));
#endif /* DO_WITH_GSL */

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


/*============================================================================*/
double randvar_normal_density(double x, double mean, double u) { 
# define CUR_PROC "randvar_normal_density"
#ifndef DO_WITH_GSL
  double expo;
#endif
  if (u <= 0.0) {
    mes_prot("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 randvar_normal_density */


/*============================================================================*/
/* special smodel 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 randvar_normal_density_approx(double x, double mean, double u) { 
# define CUR_PROC "randvar_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) {
    mes_prot("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 randvar_normal_density_approx */


/*============================================================================*/
double randvar_std_normal(int seed) {
# define CUR_PROC "randvar_std_normal"
  if (seed != 0) {
    gsl_rng_set(RNG,seed);    
    return(1.0);
  }
  else
    return(gsl_ran_gaussian(RNG, 1.0));
# undef CUR_PROC
} /* randvar_std_normal */


/*============================================================================*/
double randvar_normal(double mue, double u, int seed) {
# define CUR_PROC "randvar_normal"
#ifdef DO_WITH_GSL
  if (seed != 0) {
    gsl_rng_set(RNG,seed);    
    return(1.0*sqrt(u)+mue);
  }
  else
    return(gsl_ran_gaussian(RNG,sqrt(u))+mue);
#else
  double x;
  x = sqrt(u) * randvar_std_normal(seed) + mue;
  return(x);
#endif
# undef CUR_PROC
} /* randvar_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 randvar_normal_pos(double mue, double u, int seed) {
# define CUR_PROC "randvar_normal_pos"
  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) {
    mes_prot("u <= 0.0 not allowed\n"); goto STOP;
  }
  sigma=sqrt(u);

  if (seed != 0) {
    gsl_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, -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 * randvar_std_normal(seed) + mue;
  } */
  
  /* Inverse transformation with restricted sampling by Fishman */
  U = gsl_rng_uniform(RNG);               /*gsl_ran_flat(RNG,0,1) ??? */
  Feps = randvar_get_PHI(-(EPS_NDT+mue)/sigma);
  Us = Feps + (1-Feps)*U;
  /* Numerically better: 1-Us = 1-Feps - (1-Feps)*U, therefore: 
     Feps1 = 1-Feps, Us1 = 1-Us */
  Feps1 = randvar_get_PHI((EPS_NDT+mue)/sigma);
  Us1 = Feps1 - Feps1*U;
  t = m_min(Us,Us1);
  t = sqrt(-log(t*t));
  T = sigma * (t - (C0 + t*(C1 + t*C2)) / (1 + t*(D1 + t*(D2 + t*D3))) );
  if (Us-0.5 < 0) 
    x = mue - T;
  else
    x = mue + T;  
#endif /* DO_WITH_GSL */


STOP:
  return(x);
# undef CUR_PROC
} /* randvar_normal_pos */


/*============================================================================*/
double randvar_uniform_int(int seed, int K) {
# define CUR_PROC "randvar_uniform_int"
  if (seed != 0) {
    gsl_rng_set(RNG,seed);
    return(1.0);
  }
  else {
#ifdef DO_WITH_GSL
    /* more direct solution than old version ! */
    return (double)gsl_rng_uniform_int(RNG,K);
#else
    /* why do this ? */
    double x = gsl_ran_flat(RNG, -0.5, K-0.5);
    x = m_int(x); /* m_int rundet auf Integer */
    return(x);
#endif
  }
# undef CUR_PROC
} /* randvar_uniform_int */


/*============================================================================*/
/* cumalative distribution function of N(mean, u) */
double randvar_normal_cdf(double x, double mean, double u) { 
# define CUR_PROC "randvar_normal_cdf"
  if (u <= 0.0) {
    mes_prot("u <= 0.0 not allowed\n"); goto STOP;
  }
#ifdef DO_WITH_GSL
  /* PHI(x)=erf(x/sqrt(2))/2+0.5 */
  return (gsl_sf_erf((x-mean)/sqrt(u*2.0))+1.0)/2.0;
#else
  /* The denominator is possibly < EPS??? Check that ? */
  return(randvar_get_PHI((x - mean)/sqrt(u)));
#endif /* DO_WITH_GSL */
STOP:
  return(-1.0);
# undef CUR_PROC
} /* double randvar_normal_cdf */

/*============================================================================*/
/* cumalative distribution function of -EPS_NDT-truncated N(mean, u) */
double randvar_normal_pos_cdf(double x, double mean, double u) {
# define CUR_PROC "randvar_normal_pos_cdf"
#ifndef DO_WITH_GSL
  double Fx, c;
#endif

  if (x <= 0.0) return(0.0);
  if (u <= 0.0) { mes_prot("u <= 0.0 not allowed\n"); goto STOP; }
#ifdef DO_WITH_GSL
  /*
    Function: int gsl_sf_erfc_e (double x, gsl_sf_result * result) 
    These routines compute the complementary error function
    erfc(x) = 1 - erf(x) = 2/\sqrt(\pi) \int_x^\infty \exp(-t^2). 
   */
  return 1.0+(gsl_sf_erf((x-mean)/sqrt(u*2))-1.0)/gsl_sf_erfc((-mean)/sqrt(u*2));
#else
  /*The denominator is possibly < EPS??? Check that ? */
  Fx = randvar_get_PHI((x - mean)/sqrt(u));
  c = randvar_get_1overa(-EPS_NDT, mean, u);
  return(c*(Fx-1)+1);
#endif /* DO_WITH_GSL */
 STOP:
  return(-1.0);
# undef CUR_PROC
} /* double randvar_normal_cdf */