* samtools-0.1.8-13 (r715)

[samtools.git] / bcftools / prob1.c

diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index 85c2945779f85c9cd050fd72f00cdc78f1deeacb..b27bb8e676cbb5e45e0eb0309133d0aa7a8460e1 100644 (file)
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -106,8 +106,10 @@ int bcf_p1_read_prior(bcf_p1aux_t *ma, const char *fn)
        for (k = 0; k <= ma->M; ++k) ma->phi[k] /= sum;
        for (k = 0; k <= ma->M; ++k) fprintf(stderr, " %d:%.3lg", k, ma->phi[ma->M - k]);
        fputc('\n', stderr);
-       for (sum = 0., k = 1; k <= ma->M; ++k) sum += k * ma->phi[ma->M - k];
-       fprintf(stderr, "[heterozygosity] %lf\n", (double)sum / ma->M);
+       for (sum = 0., k = 1; k < ma->M; ++k) sum += ma->phi[ma->M - k] * (2.* k * (ma->M - k) / ma->M / (ma->M - 1));
+       fprintf(stderr, "[%s] heterozygosity=%lf, ", __func__, (double)sum);
+       for (sum = 0., k = 1; k <= ma->M; ++k) sum += k * ma->phi[ma->M - k] / ma->M;
+       fprintf(stderr, "theta=%lf\n", (double)sum);
        return 0;
 }