pssm-distr.c

EM算法的改进

/*
 * $Id: pssm-distr.c 1339 2006-09-21 19:46:28Z tbailey $
 * 
 * $Log$
 * Revision 1.2  2005/10/25 21:23:44  nadya
 * add Id and Log lines for keeping track of changes
 * change C++-style comments to C-style
 *
 *
 */

#include "macros.h"
#include "pssm-distr.h"

static double *calc_pssm_pdf(
  int w,		/* width of PSSM */
  int alen,		/* length of alphabet */
  int range,		/* largest value in PSSM */
  double **pssm,	/* scaled, integer PSSM: pssm[i][j] is score for 
			   j_th letter in i_th column of motif; 
			   entries in PSSM are in range [0..range] */
  double *prob		/* 0-order Markov background model */
);

/**********************************************************************/
/*
        calc_pssm_cdf

	Calculate the 1 minus the theoretical cumulative distribution 
	function for an (integer-valued) PSSM given a 0-order Markov
	background model.  PSSM entries must be in range [0..range].

	Returns an array of 1-cdf values:
          pvalue[x] = Pr(score >= x)
	for 0 <= x <= range*w.
*/
/**********************************************************************/
extern double *calc_pssm_cdf(
  int w,		/* width of PSSM */
  int alen,		/* length of alphabet */
  int range,		/* largest value in PSSM */
  double **pssm,	/* scaled, integer PSSM: pssm[i][j] is score for 
			   j_th letter in i_th column of motif; 
			   entries in PSSM are in range [0..range] */
  double *prob		/* 0-order Markov background model */
)
{
  int i;
  int size = w*range+1;			/* size of cdf array */
  double *pdf = NULL;	 		/* probability distribution */

  pdf = calc_pssm_pdf(w, alen, range, pssm, prob);
  if (!pdf) return NULL;

  /*for (i=0; i<size; i++) fprintf(stdout, "pdf %d %17.5e\n", i, pdf[i]);*/

  /* compute 1-cdf from the pdf from the right to preserve right accuracy */ 
  for (i=size-2; i>=0; i--) {
    pdf[i] += pdf[i+1];
    pdf[i] = MIN(1.0, pdf[i]);
    /*
      if (isnan(pdf[i])) {fprintf(stderr, "cdf: i %d pdf %f\n", i, pdf[i]);
	abort();
      }
    */
  }

  /* return the cdf */
  /*for (i=0; i<size; i++) fprintf(stdout, "cdf %d %17.5e\n", i, pdf[i]);*/
  return pdf;
} /* calc_pssm_cdf */

/**********************************************************************/
/*
        calc_pdf

	Calculate the theoretical distribution function for an 
	(integer-valued) PSSM.
	Returns an array of pdf values:
          pvalue[x] = Pr(score == x)
	for 0 <= x <= range*w.
*/
/**********************************************************************/
static double *calc_pssm_pdf(
  int w,		/* width of PSSM */
  int alen,		/* length of alphabet */
  int range,		/* largest value in PSSM */
  double **pssm,	/* scaled, integer PSSM: pssm[i][j] is score for 
			   j_th letter in i_th column of motif; 
			   entries in PSSM are in range [0..range] */
  double *prob		/* 0-order Markov background model */
)
{
  int i, j, k;
  int size = w*range+1;			/* size of pdf array */
  double *pdf_old=NULL, *pdf_new=NULL;

  /* set up the two arrays to hold probability density functions */
  Resize(pdf_old, size, double);
  Resize(pdf_new, size, double);

  /* set probabilities of each new score to zero except for score 0 */
  pdf_new[0] = 1;
  for (i=1; i<size; i++) pdf_new[i] = 0;

  /* recursively compute the pdf */
  for (i=0; i < w; i++) {		/* loop over columns in motif */
    int max_score = i * range;		/* maximum possible cumulative score */
    swap(pdf_new, pdf_old, double *);	/* new column: swap old and new pdfs */
    /* zero out the new pdf; new maximum score is old max + range */
    for (k=0; k<=max_score+range; k++) pdf_new[k] = 0;	
    for (j=0; j < alen; j++) {		/* loop over letters */
      int score = (int) pssm[i][j];	/* get integer PSSM entry */
      for (k=0; k<=max_score; k++) {
	if (pdf_old[k] != 0) pdf_new[k + score] += pdf_old[k] * prob[j];
	if (k+score >= size) {
	  fprintf(stderr, 
	    "calc_pssm_pdf error: i=%d j=%d k=%d max_score=%d score=%d size=%d\n", 
	    i, j, k, max_score, score, size);
	  return NULL;
	}
      }
    }
  }

  /* clean up */
  myfree(pdf_old);

  /* return the pdf */
  return pdf_new;

} /* calc_pssm_pdf */