Compute -log of a binomial of form p^n*(1-p)^(N-n)binomial(N,n). - Java java.lang

Description

Compute -log of a binomial of form p^n*(1-p)^(N-n)binomial(N,n).

Demo Code

/*//from  ww  w.j  a  v a2 s.com
 * Copyright (c) 2014. Real Time Genomics Limited.
 *
 * Use of this source code is bound by the Real Time Genomics Limited Software Licence Agreement
 * for Academic Non-commercial Research Purposes only.
 *
 * If you did not receive a license accompanying this file, a copy must first be obtained by email
 * from support@realtimegenomics.com.  On downloading, using and/or continuing to use this source
 * code you accept the terms of that license agreement and any amendments to those terms that may
 * be made from time to time by Real Time Genomics Limited.
 */
//package com.java2s;

public class Main {
    private static final double C0 = Math.PI / 3.0;
    private static final double C1 = Math.PI * 2.0;
    private static final int EXPONENT_OFFSET = 52;
    private static final int EXPONENT_BIAS = 1023;
    private static final double LN2 = Math.log(2);
    private static final int BITS = 16;
    private static final int MASK = (1 << BITS) - 1;
    private static final double[] LOG_TABLE = new double[1 << BITS];
    /** Cache low values to get errors lower. */
    private static final double[] LOG_F = new double[31];

    /**
     * Compute -log of a binomial of form p^n*(1-p)^(N-n)binomial(N,n).
     *
     * @param p probability
     * @param nn total population
     * @param n marked subpopulation
     * @return -log of a binomial of form p^n*(1-p)^(N-n)binomial(N,n)
     */
    public static double logBinomial(final double p, final int nn,
            final int n) {
        assert p >= 0.0 && p <= 1.0;
        assert n >= 0;
        assert nn >= 0;
        final int m = nn - n;
        if (p == 0.0) {
            if (n == 0) {
                return 0.0;
            }
            throw new IllegalArgumentException(
                    "if probability is 0.0 then count must be 0. p:" + p
                            + " N:" + nn + " n:" + n);
        }
        if (p == 1.0) {
            if (m == 0) {
                return 0.0;
            }
            throw new IllegalArgumentException(
                    "if probability is 1.0 then count must be 0. p:" + p
                            + " N:" + nn + " n:" + n);
        }
        final double res = n * Math.log(p) + m * Math.log(1.0f - p)
                + logBinomial(nn, n);
        assert res <= 0;
        return -res;
    }

    /**
     * Compute log binomial(N,n).
     *
     * @param nn total count.
     * @param n subset count.
     * @return log binomial(N,n)
     */
    public static double logBinomial(final int nn, final int n) {
        assert n >= 0;
        assert nn >= n;
        if (nn <= 1 || n == 0 || n == nn) {
            return 0.0;
        }
        if (n == 1 || n == (nn - 1)) {
            return Math.log(nn);
        }
        if (n == 2 || n == (nn - 2)) {
            return Math.log(((double) nn) * (nn - 1) / 2.0);
        }
        return logFactorial(nn) - logFactorial(n) - logFactorial(nn - n);
    }

    /**
     * Compute an approximation to the natural logarithm. Assumes
     * parameter is positive and finite.
     *
     * @param x parameter
     * @return <code>ln(x)</code>
     */
    public static double log(final double x) {
        assert x >= 0 && !Double.isInfinite(x) && !Double.isNaN(x);
        if (x == 0.0) {
            return Double.NEGATIVE_INFINITY;
        }
        final long t = Double.doubleToRawLongBits(x);
        final long lg = (t >>> EXPONENT_OFFSET) - EXPONENT_BIAS;
        final int mantissa = (int) (t >> (EXPONENT_OFFSET - BITS));
        final double mlg = LOG_TABLE[mantissa & MASK];
        return mlg + lg * LN2;
    }

    /**
     * Compute <code>logFactorial(N)</code>. Worst relative error at 31 = -9.1E-8
     *
     * @param n number
     * @return log n!
     */
    public static double logFactorial(final int n) {
        assert n >= 0;
        final double res;
        if (n < LOG_F.length) {
            res = LOG_F[n];
        } else {
            // Use Gospers Formula.
            res = Math.log(C0 + n * C1) * 0.5 + (Math.log(n) - 1.0) * n;
        }
        assert res >= 0.0 && !Double.isNaN(res) && !Double.isInfinite(res);
        return res;
    }
}

• Previous
• Next

Related Tutorials

• log Gamma
• Returns the natural log of the Binomial Coefficient, "n choose k", the number of k-element subsets that can be selected from an n-element set.
• Given log(a) and log(b) calculate log(a + b) boils down to log( exp(log_a) + exp(log_b) ) but this might overflow, so we turn this into log([exp(log_a - log_c) + exp(log_b - log_c)]exp(log_c)) and we set log_c == max(log_a,log_b) and so it becomes: LARGE + log(1 + exp(SMALL - LARGE)) == LARGE + log(1 + SMALL) ~= large the whole idea being to avoid an overflow (exp(LARGE) == VERY LARGE == overflow)
• turn a log probability into Phred score
• Compute an approximation to the natural logarithm.