--- /dev/null
+/* Ethan Wolf 1995 */
+
+#include <stdio.h>
+#include <math.h>
+#include "methods.h"
+
+
+
+/*** Note that to get the gaussian on the same scale as the histogram of **/
+/*** residue counts, we need to take counts[f][i]/number_residues_in_hist. **/
+double gauss_value(double x, double mean, double sd, double down)
+{
+ double gauss_x;
+
+ /*** Note: 2*pi=2.506... **/
+
+ gauss_x = down * exp( - (x-mean)/sd*(x-mean)/sd /2 ) /
+ (sd*2.5066282746310002);
+ return(gauss_x);
+}
+
+/** The following approximates the gaussian by the line from gauss1 to */
+/** gauss2. Let y= hist(x). y1=min{gauss1,gauss2}, y2=max{gauss1,gauss2} */
+/** Case1:
+ y>y2
+ Then area of hist over gauss is:
+ y-y2 + .5*(y2-y1)
+ Case2:
+ y<y1
+ THen area of hist over gauss is 0.
+ Case3:
+ y1<y<y2
+ Then fraction width fraction (y-y1)/(y2-y1) is over gauss.
+ so that contributes .5 * (y-y1)/(y2-y1) * (y-y1) to area.
+**********************/
+
+
+
+
+
+
+double like_gauss_by_line(double x, double mean, double sd, double down,
+ double hist_height)
+{
+ double gauss1, gauss2, y1, y2;
+ double over_gauss;
+
+ /*** Just put this in on Jan 15, 1996 cuz I don't think we want to use */
+ /*** the down value, since the histogram was just fit to the left */
+ /*** side of scores below the mean. *************/
+
+/* down=1; */
+ /*******************/
+
+ gauss1 = gauss_value(x-.25,mean,sd,down);
+ gauss2= gauss_value(x+.25,mean,sd,down);
+
+ if (gauss1<gauss2) {
+ y1=gauss1; y2=gauss2; }
+ else {
+ y2=gauss1; y1=gauss2; }
+
+ if (hist_height>=y2)
+ return( (hist_height-y2 + .5*(y2-y1))/hist_height);
+ else if (hist_height<=y1)
+ return(0);
+ else {
+ over_gauss = .5* ( (hist_height-y1)*(hist_height-y1)/(y2-y1) -
+ - (y2- hist_height)*(y2- hist_height)/(y2-y1));
+ if (over_gauss>0)
+ return(over_gauss/hist_height);
+ else
+ return(0);
+ }
+}
+
+
+
+/*** The following considers the area under the gauss to be
+ gauss(x) so the height of the gaussian at x (x could be
+ the midpoint of the box, or the mean of the scores
+ in the box, etc.
+****/
+
+double like_gauss_by_box(double x, double mean, double sd, double down,
+ double hist_height)
+{
+ double gauss;
+
+
+ gauss = gauss_value(x,mean,sd,down);
+ if (hist_height >= gauss)
+ return( (hist_height - gauss)/hist_height);
+ else
+ return(0);
+}
+
+
+
+/*** The following function takes two gaussians (pdb- and the positives) **/
+/*** and a parameter (estimate of the ratio of negatives to positives in **/
+/*** all proteins) and computes likelihood(x) as:
+ gauss_pos(x)/[gauss_pos(x) + ratio*gauss_neg(x)]
+**/
+/* This is the method that STOCK used for likelihoods. **/
+
+double stock_like(double ratio, double x, double mean_neg, double sd_neg,
+ double down_neg, double mean_pos, double sd_pos, double down_pos)
+{
+ double gauss_neg, gauss_pos;
+
+ gauss_neg= gauss_value(x,mean_neg,sd_neg,down_neg);
+ gauss_pos= gauss_value(x,mean_pos,sd_pos,down_neg);
+
+ return(gauss_pos/(gauss_pos + ratio*gauss_neg));
+}
+
+
+/**** The following function computes likelihood(x) as:
+ fraction of residues that score AT LEAST x that are "positives"
+ i.e area_pos_score_at_least_x/area_total_score_at_least_x
+
+ Use total_area_above_gauss() and total_area() to compute initial
+ values for area_pos_score_at_least_x and area_all_score_at_least_x.
+****/
+
+double ratio_PosToNeg_at_least_x(double x, double mean, double sd, double down,
+ double *area_pos_score_at_least_x,
+ double *area_all_score_at_least_x,
+ double hist_height)
+{
+ double gauss, like;
+
+ like = (*area_pos_score_at_least_x)/(*area_all_score_at_least_x);
+
+ gauss = gauss_value(x,mean,sd,down);
+ if (hist_height > gauss)
+ *area_pos_score_at_least_x -= hist_height -gauss;
+
+ *area_all_score_at_least_x -= hist_height;
+
+ return(like);
+}
+
+
+/*** This function is not so much a likelihood, as the fraction of
+ "positives" (things above the gaussian) that score x or less. ***/
+/*** To compute this we need to first run through and compute the total **/
+/*** area above the gaussian. We also keep a running total of the area **/
+/*** above the gaussian up to x. **/
+
+double fraction_of_positives_lower_than_x(double x, double mean, double sd,
+ double down, double *area_to_x, double total_area,
+ double hist_height)
+{
+ double gauss;
+
+ gauss = gauss_value(x,mean,sd,down);
+ if (hist_height > gauss)
+ *area_to_x += hist_height - gauss;
+
+ return(*area_to_x/total_area);
+}
+
+
+
+
+
+/*** hist_height[i] corresponds to the x value i/2 ***/
+double total_area_above_gauss(double min_x, double max_x, long *counts,
+ double mean, double sd, double down,
+ long number_residues)
+{
+ double hist_height;
+ double gauss, area_above=0;
+ int i;
+
+ for(i=2*min_x; i<= 2*max_x; i++) {
+ gauss = gauss_value((double)i/2,mean,sd,down);
+ /** Need to convert counts to be on same scale as gauss. **/
+ hist_height = 2*(double)counts[i]/number_residues;
+
+ if (hist_height > gauss)
+ area_above += hist_height-gauss;
+ }
+
+ return(area_above);
+}
+
+double total_area(double min_x, double max_x, long *counts,
+ long number_residues)
+{
+ double hist_height;
+ double area=0;
+ int i;
+
+ for(i=2*min_x; i<= 2*max_x; i++) {
+ /** Need to convert counts to be on same scale as gauss. **/
+ hist_height = 2*(double)counts[i]/number_residues;
+ area += hist_height;
+ }
+ return(area);
+}
git://source.jalview.org / jpred.git / blobdiff