#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
//____________________________________________________________________________