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)
15 /*** Note: 2*pi=2.506... **/
17 gauss_x = down * exp( - (x-mean)/sd*(x-mean)/sd /2 ) /
18 (sd*2.5066282746310002);
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} */
26 Then area of hist over gauss is:
30 THen area of hist over gauss is 0.
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 **********************/
42 double like_gauss_by_line(double x, double mean, double sd, double down,
45 double gauss1, gauss2, y1, y2;
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. *************/
55 gauss1 = gauss_value(x-.25,mean,sd,down);
56 gauss2= gauss_value(x+.25,mean,sd,down);
59 y1=gauss1; y2=gauss2; }
61 y2=gauss1; y1=gauss2; }
64 return( (hist_height-y2 + .5*(y2-y1))/hist_height);
65 else if (hist_height<=y1)
68 over_gauss = .5* ( (hist_height-y1)*(hist_height-y1)/(y2-y1) -
69 - (y2- hist_height)*(y2- hist_height)/(y2-y1));
71 return(over_gauss/hist_height);
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
85 double like_gauss_by_box(double x, double mean, double sd, double down,
91 gauss = gauss_value(x,mean,sd,down);
92 if (hist_height >= gauss)
93 return( (hist_height - gauss)/hist_height);
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)]
105 /* This is the method that STOCK used for likelihoods. **/
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)
110 double gauss_neg, gauss_pos;
112 gauss_neg= gauss_value(x,mean_neg,sd_neg,down_neg);
113 gauss_pos= gauss_value(x,mean_pos,sd_pos,down_neg);
115 return(gauss_pos/(gauss_pos + ratio*gauss_neg));
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
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.
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,
134 like = (*area_pos_score_at_least_x)/(*area_all_score_at_least_x);
136 gauss = gauss_value(x,mean,sd,down);
137 if (hist_height > gauss)
138 *area_pos_score_at_least_x -= hist_height -gauss;
140 *area_all_score_at_least_x -= hist_height;
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. **/
152 double fraction_of_positives_lower_than_x(double x, double mean, double sd,
153 double down, double *area_to_x, double total_area,
158 gauss = gauss_value(x,mean,sd,down);
159 if (hist_height > gauss)
160 *area_to_x += hist_height - gauss;
162 return(*area_to_x/total_area);
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)
175 double gauss, area_above=0;
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;
183 if (hist_height > gauss)
184 area_above += hist_height-gauss;
190 double total_area(double min_x, double max_x, long *counts,
191 long number_residues)
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;