sources/multicoil/methods.c

   1 /* Ethan Wolf 1995 */
   2
   3 #include <stdio.h>
   4 #include <math.h>
   5 #include "methods.h"
   6
   7
   8
   9 /*** Note that to get the gaussian on the same scale as the histogram of   **/
  10 /*** residue counts, we need to take counts[f][i]/number_residues_in_hist. **/
  11 double gauss_value(double x, double mean, double sd, double down)
  12 {
  13   double gauss_x;
  14
  15    /*** Note: 2*pi=2.506... **/
  16
  17   gauss_x = down * exp( - (x-mean)/sd*(x-mean)/sd /2 ) /
  18                   (sd*2.5066282746310002);
  19   return(gauss_x);
  20 }
  21
  22 /** The following approximates the gaussian by the line from gauss1 to */
  23 /** gauss2.  Let y= hist(x). y1=min{gauss1,gauss2}, y2=max{gauss1,gauss2}    */
  24 /** Case1:
  25        y>y2
  26            Then area of hist over gauss is:
  27            y-y2 + .5*(y2-y1)
  28     Case2:
  29        y<y1
  30            THen area of hist over gauss is 0.
  31     Case3:
  32        y1<y<y2
  33            Then fraction width fraction (y-y1)/(y2-y1) is over  gauss.
  34              so that contributes  .5 * (y-y1)/(y2-y1) *  (y-y1)  to area.
  35 **********************/
  36
  37
  38
  39
  40
  41
  42 double like_gauss_by_line(double x, double mean, double sd, double down,
  43                      double hist_height)
  44 {
  45   double gauss1, gauss2, y1, y2;
  46   double over_gauss;
  47
  48   /*** Just put this in on Jan 15, 1996 cuz I don't think we want to use */
  49   /*** the down value, since the histogram was just fit to the left */
  50   /*** side of scores below the mean.  *************/
  51
  52 /*  down=1; */
  53  /*******************/
  54
  55   gauss1 = gauss_value(x-.25,mean,sd,down);
  56   gauss2= gauss_value(x+.25,mean,sd,down);
  57
  58   if (gauss1<gauss2) {
  59     y1=gauss1;  y2=gauss2; }
  60   else {
  61     y2=gauss1; y1=gauss2;  }
  62
  63   if (hist_height>=y2)
  64     return( (hist_height-y2 + .5*(y2-y1))/hist_height);
  65   else if (hist_height<=y1)
  66     return(0);
  67   else {
  68     over_gauss = .5* ( (hist_height-y1)*(hist_height-y1)/(y2-y1)  -
  69                       - (y2- hist_height)*(y2- hist_height)/(y2-y1));
  70     if (over_gauss>0)
  71       return(over_gauss/hist_height);
  72     else
  73       return(0);
  74   }
  75 }
  76
  77
  78
  79 /*** The following considers the area under the gauss to be
  80      gauss(x) so the height of the gaussian at x (x could be
  81      the midpoint of the box, or the mean of the scores
  82      in the box, etc.
  83 ****/
  84
  85 double like_gauss_by_box(double x, double mean, double sd, double down,
  86                          double hist_height)
  87 {
  88   double gauss;
  89
  90
  91   gauss = gauss_value(x,mean,sd,down);
  92   if (hist_height >= gauss)
  93     return( (hist_height - gauss)/hist_height);
  94   else
  95     return(0);
  96 }
  97
  98
  99
 100 /*** The following function takes two gaussians (pdb- and the positives) **/
 101 /*** and a parameter (estimate of the ratio of negatives to positives in **/
 102 /*** all proteins) and computes likelihood(x) as:
 103           gauss_pos(x)/[gauss_pos(x) + ratio*gauss_neg(x)]
 104 **/
 105 /*  This is the method that STOCK used for likelihoods.                 **/
 106
 107 double stock_like(double ratio, double x, double mean_neg, double sd_neg,
 108              double down_neg, double mean_pos, double sd_pos, double down_pos)
 109 {
 110   double gauss_neg, gauss_pos;
 111
 112   gauss_neg=  gauss_value(x,mean_neg,sd_neg,down_neg);
 113   gauss_pos=  gauss_value(x,mean_pos,sd_pos,down_neg);
 114
 115   return(gauss_pos/(gauss_pos + ratio*gauss_neg));
 116 }
 117
 118
 119 /**** The following function computes likelihood(x) as:
 120   fraction of residues that score AT LEAST x that are "positives"
 121  i.e area_pos_score_at_least_x/area_total_score_at_least_x
 122
 123  Use total_area_above_gauss() and total_area() to compute initial
 124  values for area_pos_score_at_least_x and area_all_score_at_least_x.
 125 ****/
 126
 127 double ratio_PosToNeg_at_least_x(double x, double mean, double sd, double down,
 128                                double *area_pos_score_at_least_x,
 129                                double *area_all_score_at_least_x,
 130                                double hist_height)
 131 {
 132   double gauss, like;
 133
 134   like = (*area_pos_score_at_least_x)/(*area_all_score_at_least_x);
 135
 136   gauss = gauss_value(x,mean,sd,down);
 137   if (hist_height > gauss)
 138     *area_pos_score_at_least_x -= hist_height -gauss;
 139
 140   *area_all_score_at_least_x -= hist_height;
 141
 142   return(like);
 143 }
 144
 145
 146 /*** This function is not so much a likelihood, as the fraction of
 147      "positives" (things above the gaussian) that score x or less. ***/
 148 /***  To compute this we need to first run through and compute the total **/
 149 /***  area above the gaussian.  We also keep a running total of the area **/
 150 /***  above the gaussian up to x.   **/
 151
 152 double fraction_of_positives_lower_than_x(double x, double mean, double sd,
 153                   double down, double *area_to_x, double total_area,
 154                   double hist_height)
 155 {
 156   double gauss;
 157
 158   gauss = gauss_value(x,mean,sd,down);
 159   if (hist_height > gauss)
 160     *area_to_x += hist_height - gauss;
 161
 162   return(*area_to_x/total_area);
 163 }
 164
 165
 166
 167
 168
 169 /*** hist_height[i] corresponds to the x value i/2 ***/
 170 double total_area_above_gauss(double min_x, double max_x, long *counts,
 171                               double mean, double sd, double down,
 172                               long number_residues)
 173 {
 174   double hist_height;
 175   double gauss, area_above=0;
 176   int i;
 177
 178   for(i=2*min_x; i<= 2*max_x; i++) {
 179     gauss = gauss_value((double)i/2,mean,sd,down);
 180     /** Need to convert counts to be on same scale as gauss.   **/
 181     hist_height = 2*(double)counts[i]/number_residues;
 182
 183     if (hist_height > gauss)
 184       area_above += hist_height-gauss;
 185   }
 186
 187   return(area_above);
 188 }
 189
 190 double total_area(double min_x, double max_x, long *counts,
 191                   long number_residues)
 192 {
 193   double hist_height;
 194   double area=0;
 195   int i;
 196
 197   for(i=2*min_x; i<= 2*max_x; i++) {
 198     /** Need to convert counts to be on same scale as gauss.   **/
 199     hist_height = 2*(double)counts[i]/number_residues;
 200     area += hist_height;
 201   }
 202   return(area);
 203 }