minaux1_t *a = (minaux1_t*)data;
double p = 1., l = 0., f3[3];
int i;
-// printf("%lg\n", f);
+// printf("brent %lg\n", f);
if (f < 0 || f > 1) return 1e300;
f3[0] = (1.-f)*(1.-f); f3[1] = 2.*f*(1.-f); f3[2] = f*f;
for (i = a->beg; i < a->end; ++i) {
{
double f0 = *f, f3[3], err;
int i;
+// printf("em %lg\n", *f);
f3[0] = (1.-f0)*(1.-f0); f3[1] = 2.*f0*(1.-f0); f3[2] = f0*f0;
for (i = beg, f0 = 0.; i < end; ++i) {
const double *pdg = _pdg + i * 3;
// x[5..6]: group1 freq, group2 freq
// x[7]: 1-degree P-value
// x[8]: 2-degree P-value
-int bcf_em1(const bcf1_t *b, int n1, int flag, double x[9])
+// x[9]: 1-degree P-value with freq estimated from genotypes
+int bcf_em1(const bcf1_t *b, int n1, int flag, double x[10])
{
double *pdg;
int i, n, n2;
if (flag & 0xf<<1) flag |= 0xf<<1;
n = b->n_smpl; n2 = n - n1;
pdg = get_pdg3(b);
- for (i = 0; i < 9; ++i) x[i] = -1.;
+ if (pdg == 0) return -1;
+ for (i = 0; i < 10; ++i) x[i] = -1.; // set to negative
{
if ((x[0] = est_freq(n, pdg)) < 0.) {
free(pdg);
}
x[0] = freqml(x[0], 0, n, pdg);
}
- if (flag & (0xf<<1|1<<8)) { // estimate the genotype frequency and test HWE
+ if (flag & (0xf<<1|3<<8)) { // estimate the genotype frequency and test HWE
double *g = x + 1, f3[3], r;
f3[0] = g[0] = (1 - x[0]) * (1 - x[0]);
f3[1] = g[1] = 2 * x[0] * (1 - x[0]);
f3[i][0] = (1-f[i])*(1-f[i]), f3[i][1] = 2*f[i]*(1-f[i]), f3[i][2] = f[i]*f[i];
x[7] = kf_gammaq(.5, log(lk_ratio_test(n, n1, pdg, f3)));
}
- if ((flag & 1<<8) && n1 > 0 && n1 < n) { // 2-degree P-value
- double g[3][3];
+ if ((flag & 3<<8) && n1 > 0 && n1 < n) { // 2-degree P-value
+ double g[3][3], tmp;
for (i = 0; i < 3; ++i) memcpy(g[i], x + 1, 3 * sizeof(double));
for (i = 0; i < ITER_MAX; ++i)
if (g3_iter(g[1], pdg, 0, n1) < EPS) break;
for (i = 0; i < ITER_MAX; ++i)
if (g3_iter(g[2], pdg, n1, n) < EPS) break;
- x[8] = kf_gammaq(1., log(lk_ratio_test(n, n1, pdg, g)));
+ tmp = log(lk_ratio_test(n, n1, pdg, g));
+ x[8] = kf_gammaq(1., tmp);
+ x[9] = kf_gammaq(.5, tmp);
}
// free
free(pdg);