#include "includes.h" #include "grid.h" #include "vector.h" /* This file contains source code for functions that generate the confidence intervals for the location of the gene */ //____________________________________________________________________________ // C Routine to calculate Mary Sara's smapprox function for a specified n double smapprox(int n) { double* LGAMMA = new double[n + 100001]; LGAMMA[0] = LGAMMA[1] = 0; for (int i = 2; i < n + 100001; ++i) LGAMMA[i] = LGAMMA[i - 1] + log(i - 1); double res = 0; for (int k = 1; k <= n - 2; ++k) { double need = 0; for (int c = 0; c <= 100000; ++c) //was c = 1 need += pow(-1.0, c) * exp(LGAMMA[k + c] - LGAMMA[k] + LGAMMA[n] - LGAMMA[n + c]); res += exp(log(2) + log(n + 1) + log(n - k - 1) - 2 * log(n - 1) - log(n - k + 1) - log(n - k + 2)) * need; } delete [] LGAMMA; return res; } // smapprox //____________________________________________________________________________ // Given 'npts', the number of locations, 'd', the vector of locations, // 'll', the vector of log-likelihoods at each location, and 'chip', a // chi-square (1 df) percentile, return a list of intervals containing // the "likely" true locations. The values returned are 'edge' and // 'numedge'. 'numedge' is the number of intervals, and 'edge' is a // vector of length 2 * 'numedge' that contains pairs of left and right // indices in 'd'. The following code prints out the edge intervals. // for (int i = 0; i < numedge; ++i) // cout << "[" << d[edge[2 * i]] << ", " << d[edge[2 * i + 1]] << "]\n"; void cr(int npts, const dvector& d, const dvector& ll, double chip, ivector& edge, int& numedge) { /* KJW - If the code to calculate chi-square quantiles is not desired, we could use this array of common quantiles instead. I found the code to calculate the quantiles on netlib, modified it, and I believe there is no problem using it here. static double chipctl[12] = {.75, .80, .85, .90, .95, .975, .98, .99, .995, .9975, .999, .9995 }; static double chiqntl[12] = {1.32, 1.64, 2.07, 2.71, 3.84, 5.02, 5.41, 6.63, 7.88, 9.14, 10.83, 12.12 }; for (int i = 0; i < 12; ++i) if (chip == chipctl[i]) chiq = chiqntl[i]; */ double critchi (double p); double maxll = ll[ll.maxind()]; double chiq = critchi(1 - chip); double th = maxll - chiq / 2.0; ivector maxedge(1000); int currpt = 0; int curredge = 0; while (currpt < npts) { while (currpt < npts && ll[currpt] < th) ++currpt; if (currpt < npts) { maxedge[curredge++] = currpt; while (currpt < npts && ll[currpt] >= th) ++currpt; maxedge[curredge++] = currpt - 1; } } numedge = curredge / 2; edge.init(curredge); for (int i = 0; i < curredge; ++i) edge[i] = maxedge[i]; } // cr //____________________________________________________________________________ /* Code modified by Kenneth Wilder. Relevant original credits follow here and elsewhere in the code. This code is only used in the 'critchi' call in the 'cr' function. Programmer: Gary Perlman Organization: Wang Institute, Tyngsboro, MA 01879 Copyright: none */ #define CHI_EPSILON 0.000001 /* accuracy of critchi approximation */ #define CHI_MAX 99999.0 /* maximum chi square value */ #define LOG_SQRT_PI 0.5723649429247000870717135 /* log (sqrt (pi)) */ #define I_SQRT_PI 0.5641895835477562869480795 /* 1 / sqrt (pi) */ #define BIGX 20.0 /* max value to represent exp (x) */ #define ex(x) (((x) < -BIGX) ? 0.0 : exp (x)) double pochisq (double x); double critchi (double p); double poz (double z); //____________________________________________________________________________ // Return x so that P(X>x)='p', X a chi-square RV with 1 df. // Uses a simple bisection method. double critchi (double p) { double minchisq = 0.0; double maxchisq = CHI_MAX; if (p <= 0.0) return (maxchisq); else if (p >= 1.0) return (0.0); double chisqval = 1.0 / sqrt (p); /* fair first value */ while (maxchisq - minchisq > CHI_EPSILON) { if (pochisq(chisqval) < p) maxchisq = chisqval; else minchisq = chisqval; chisqval = (maxchisq + minchisq) * 0.5; } return (chisqval); } // critchi //____________________________________________________________________________ // Return P(X>'x')=p, X a chi-square RV with 1 df. /*FUNCTION pochisq: probability of chi sqaure value Adapted from: Hill, I. D. and Pike, M. C. Algorithm 299 Collected Algorithms for the CACM 1967 p. 243 Updated for rounding errors based on remark in ACM TOMS June 1985, page 185 */ double pochisq (double x) { if (x <= 0.0) return (1.0); else return 2.0 * poz(-sqrt(x)); } // pochisq //____________________________________________________________________________ // Return P(Z<'z'), Z a N(0,1) RV. /*FUNCTION poz: probability of normal z value Adapted from a polynomial approximation in: Ibbetson D, Algorithm 209 Collected Algorithms of the CACM 1963 p. 616 Note: This routine has six digit accuracy, so it is only useful for absolute z values < 6. For z values >= to 6.0, poz() returns 0.0. */ double poz (double z) { double x = 0.0; if (z != 0.0) { double y = 0.5 * fabs (z); if (y >= 3.0) x = 1.0; else if (y < 1.0) { double w = y*y; x = ((((((((0.000124818987 * w -0.001075204047) * w +0.005198775019) * w -0.019198292004) * w +0.059054035642) * w -0.151968751364) * w +0.319152932694) * w -0.531923007300) * w +0.797884560593) * y * 2.0; } else { y -= 2.0; x = (((((((((((((-0.000045255659 * y +0.000152529290) * y -0.000019538132) * y -0.000676904986) * y +0.001390604284) * y -0.000794620820) * y -0.002034254874) * y +0.006549791214) * y -0.010557625006) * y +0.011630447319) * y -0.009279453341) * y +0.005353579108) * y -0.002141268741) * y +0.000535310849) * y +0.999936657524; } } return (z > 0.0 ? ((x + 1.0) * 0.5) : ((1.0 - x) * 0.5)); } // poz //____________________________________________________________________________