Probability Calculations Numerically Stable

ProbabilityUtils.java

package processing.utils;

/**
 *   Utilities to compute log probabilities
 *
 *   @author Daniel Kohlsdorf
 */
public class ProbabilityUtils {
    
    public static final double ZERO = Double.NEGATIVE_INFINITY;

    public static double SCALER = Math.log(1.0 / Math.sqrt(2 * Math.PI));
    
    /**
     * Log of a data point, ZERO if x = 0
     */
    public static double lg(double x) {
        if(x == 0) {
            return ZERO;
        }
        return Math.log(x);
    }

    /**
     * Multivariate spherical normal distribution
     *
     * @param x multi dimensional point
     * @param mean multi dimensional mean
     * @param variance scalar variance
     *
     * @return log N(x | mu, sigma)
     */
    public static double lgnormpdf(double x[], double[] mean, double variance[]) {
        double ll = 0;
        for(int i = 0; i < x.length; i++) {
            double pdf = SCALER - 0.5 * lg(variance[i]) - (Math.pow(x[i] - mean[i], 2) / (2.0 * variance[i]));
            ll += pdf;
        }
        return ll;
    }
    
    /**
     * Numerically stable sum(log(x), log(y))
     *
     * @param x
     * @param y
     * @return sum(log(x), log(y))
     */
    public static double sum(double log_x, double log_y) {
        if (log_x == ZERO || log_y == ZERO) {
            if (log_x == ZERO) {
                // y + 0 = y
                return log_y;
            } else {
                // x + 0 = x
                return log_x;
            }
        } else {
            if (log_x > log_y) {
                return log_x + Math.log(1 + Math.exp(log_y - log_x));
            } else {
                return log_y + Math.log(1 + Math.exp(log_x - log_y));
            }
        }
    }

    /**
     * Numerical stable exp
     *
     * @param log_x
     * @return exp(log_x)
     */
    public static double exp(double log_x) {
        if (log_x == ZERO) {
            return 0;
        }
        return Math.exp(log_x);
    }

    
}