};
struct __bcf_p1aux_t {
- int n, M, n1, is_indel;
+ int n, M, n1, is_indel, is_folded;
double *q2p, *pdg; // pdg -> P(D|g)
double *phi, *phi_indel;
double *z, *zswap; // aux for afs
int PL_len;
};
+static void fold_array(int M, double *x)
+{
+ int k;
+ for (k = 0; k < M/2; ++k)
+ x[k] = x[M-k] = (x[k] + x[M-k]) / 2.;
+}
+
void bcf_p1_indel_prior(bcf_p1aux_t *ma, double x)
{
int i;
return 0;
}
+void bcf_p1_set_folded(bcf_p1aux_t *p1a)
+{
+ if (p1a->n1 < 0) {
+ p1a->is_folded = 1;
+ fold_array(p1a->M, p1a->phi);
+ fold_array(p1a->M, p1a->phi_indel);
+ }
+}
+
void bcf_p1_destroy(bcf_p1aux_t *ma)
{
if (ma) {
//
rst->rank0 = cal_pdg(b, ma);
rst->f_exp = mc_cal_afs(ma);
- rst->p_ref = ma->afs1[ma->M];
+ rst->p_ref = ma->is_folded? ma->afs1[ma->M] + ma->afs1[0] : ma->afs1[ma->M];
// calculate f_flat and f_em
for (k = 0, sum = 0.; k <= ma->M; ++k)
sum += (long double)ma->z[k];
flast = rst->f_em;
}
}
+ { // estimate equal-tail credible interval (95% level)
+ int l, h;
+ double p;
+ for (i = 0, p = 0.; i < ma->M; ++i)
+ if (p + ma->afs1[i] > 0.025) break;
+ else p += ma->afs1[i];
+ l = i;
+ for (i = ma->M-1, p = 0.; i >= 0; --i)
+ if (p + ma->afs1[i] > 0.025) break;
+ else p += ma->afs1[i];
+ h = i;
+ rst->cil = (double)(ma->M - h) / ma->M; rst->cih = (double)(ma->M - l) / ma->M;
+ }
rst->g[0] = rst->g[1] = rst->g[2] = -1.;
contrast(ma, rst->pc);
return 0;
void bcf_p1_dump_afs(bcf_p1aux_t *ma)
{
int k;
+ if (ma->is_folded) fold_array(ma->M, ma->afs);
fprintf(stderr, "[afs]");
for (k = 0; k <= ma->M; ++k)
fprintf(stderr, " %d:%.3lf", k, ma->afs[ma->M - k]);