www.pudn.com > tpWY.rar > stat_util.c

/* Miscellaneous utility statistical functions. */

/* Created: Tue Oct 24 15:54:19 EDT 2000 */
/* Elisabetta Manduchi, CBIL, Penn Center for Bioinformatics */

#include 
#include 
#include "stat_util.h"

/* Returns the binomial coefficient "n chooose k" as a floating-point 
number. */
float bicoeff(int n, int k) {
  return floor(0.5+exp(lnfact(n)-lnfact(k)-lnfact(n-k)));
}

/* Returns ln(n!). */
float lnfact(int n) {
  if (n<0) {
    fprintf(stderr, "Negative factorial in routine lnfact.\n");
    exit(1);
  }
  if (n<=1)
    return 0.0;
  else
    return lngamma(n+1.0);
}

/* Returns the value ln|Gamma(z)| for z>0 (uses Lanczos' approximation). */
float lngamma(float z) {
  double x, lng, series;
  static double coeff[6] = {76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5};
  int i;

  x = z;
  lng = -z-5.5;
  lng += (z+0.5)*log(-lng);
  series = 1.000000000190015;
  for (i=0; i<5; i++)
    series += coeff[i]/++x;
  return lng+log(2.5066282746310005*series/z);
}

/* Given a list a of k numbers a_0=0)
    *(ptr+h) += 1;
}

/* Returns a random number between 0 and 1 */
float ran2(long *idum)
{
  static long iy,ir[98];
  static int iff=0;
  int j;
  
  if (*idum < 0 || iff == 0) {
    iff=1;
    if ((*idum=(IC-(*idum)) % M) < 0) *idum = -(*idum);
    for (j=1; j<=97; j++) {
      *idum=(IA*(*idum)+IC) % M;
      ir[j]=(*idum);
    }
    *idum=(IA*(*idum)+IC) % M;
    iy=(*idum);
  }
  j=1 + 97.0*iy/M;
  if (j > 97 || j < 1) {
    fprintf(stderr, "RAN2: This cannot happen.\n");
  }
  iy=ir[j];
  *idum=(IA*(*idum)+IC) % M;
  ir[j]=(*idum);
  return (float) iy/M;
}

/* Generates a random subset of k elements from the set {0, 1,..., n-1} and
places it into *ptr */
void ran_subset(int *ptr, int n, int k, long *s) {
  int *arr, i, j;

  arr = (int *) malloc(n * sizeof(int)); 
  if (arr==NULL) {
    fprintf(stderr, "insufficient memory (ran_subset::arr)\n");
  }
  for (i=0; i