FoldedNormalDistribution.java

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.commons.statistics.distribution;

import org.apache.commons.numbers.gamma.Erf;
import org.apache.commons.numbers.gamma.ErfDifference;
import org.apache.commons.numbers.gamma.Erfc;
import org.apache.commons.numbers.gamma.InverseErf;
import org.apache.commons.numbers.gamma.InverseErfc;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.sampling.distribution.GaussianSampler;
import org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler;
import org.apache.commons.rng.sampling.distribution.ZigguratSampler;

/**
 * Implementation of the folded normal distribution.
 *
 * <p>Given a normally distributed random variable \( X \) with mean \( \mu \) and variance
 * \( \sigma^2 \), the random variable \( Y = |X| \) has a folded normal distribution. This is
 * equivalent to not recording the sign from a normally distributed random variable.
 *
 * <p>The probability density function of \( X \) is:
 *
 * <p>\[ f(x; \mu, \sigma) = \frac 1 {\sigma\sqrt{2\pi}} e^{-{\frac 1 2}\left( \frac{x-\mu}{\sigma} \right)^2 } +
 *                           \frac 1 {\sigma\sqrt{2\pi}} e^{-{\frac 1 2}\left( \frac{x+\mu}{\sigma} \right)^2 }\]
 *
 * <p>for \( \mu \) the location,
 * \( \sigma &gt; 0 \) the scale, and
 * \( x \in [0, \infty) \).
 *
 * <p>If the location \( \mu \) is 0 this reduces to the half-normal distribution.
 *
 * @see <a href="https://en.wikipedia.org/wiki/Folded_normal_distribution">Folded normal distribution (Wikipedia)</a>
 * @see <a href="https://en.wikipedia.org/wiki/Half-normal_distribution">Half-normal distribution (Wikipedia)</a>
 * @since 1.1
 */
public abstract class FoldedNormalDistribution extends AbstractContinuousDistribution {
    /** The scale. */
    final double sigma;
    /**
     * The scale multiplied by sqrt(2).
     * This is used to avoid a double division when computing the value passed to the
     * error function:
     * <pre>
     *  ((x - u) / sd) / sqrt(2) == (x - u) / (sd * sqrt(2)).
     *  </pre>
     * <p>Note: Implementations may first normalise x and then divide by sqrt(2) resulting
     * in differences due to rounding error that show increasingly large relative
     * differences as the error function computes close to 0 in the extreme tail.
     */
    final double sigmaSqrt2;
    /**
     * The scale multiplied by sqrt(2 pi). Computed to high precision.
     */
    final double sigmaSqrt2pi;

    /**
     * Regular implementation of the folded normal distribution.
     */
    private static class RegularFoldedNormalDistribution extends FoldedNormalDistribution {
        /** The location. */
        private final double mu;
        /** Cached value for inverse probability function. */
        private final double mean;
        /** Cached value for inverse probability function. */
        private final double variance;

        /**
         * @param mu Location parameter.
         * @param sigma Scale parameter.
         */
        RegularFoldedNormalDistribution(double mu, double sigma) {
            super(sigma);
            this.mu = mu;

            final double a = mu / sigmaSqrt2;
            mean = sigma * Constants.ROOT_TWO_DIV_PI * Math.exp(-a * a) + mu * Erf.value(a);
            this.variance = mu * mu + sigma * sigma - mean * mean;
        }

        @Override
        public double getMu() {
            return mu;
        }

        @Override
        public double density(double x) {
            if (x < 0) {
                return 0;
            }
            final double vm = (x - mu) / sigma;
            final double vp = (x + mu) / sigma;
            return (ExtendedPrecision.expmhxx(vm) + ExtendedPrecision.expmhxx(vp)) / sigmaSqrt2pi;
        }

        @Override
        public double probability(double x0,
                                  double x1) {
            if (x0 > x1) {
                throw new DistributionException(DistributionException.INVALID_RANGE_LOW_GT_HIGH,
                                                x0, x1);
            }
            if (x0 <= 0) {
                return cumulativeProbability(x1);
            }
            // Assumes x1 >= x0 && x0 > 0
            final double v0m = (x0 - mu) / sigmaSqrt2;
            final double v1m = (x1 - mu) / sigmaSqrt2;
            final double v0p = (x0 + mu) / sigmaSqrt2;
            final double v1p = (x1 + mu) / sigmaSqrt2;
            return 0.5 * (ErfDifference.value(v0m, v1m) + ErfDifference.value(v0p, v1p));
        }

        @Override
        public double cumulativeProbability(double x) {
            if (x <= 0) {
                return 0;
            }
            return 0.5 * (Erf.value((x - mu) / sigmaSqrt2) + Erf.value((x + mu) / sigmaSqrt2));
        }

        @Override
        public double survivalProbability(double x) {
            if (x <= 0) {
                return 1;
            }
            return 0.5 * (Erfc.value((x - mu) / sigmaSqrt2) + Erfc.value((x + mu) / sigmaSqrt2));
        }

        @Override
        public double getMean() {
            return mean;
        }

        @Override
        public double getVariance() {
            return variance;
        }

        @Override
        public Sampler createSampler(UniformRandomProvider rng) {
            // Return the absolute of a Gaussian distribution sampler.
            final SharedStateContinuousSampler s =
                GaussianSampler.of(ZigguratSampler.NormalizedGaussian.of(rng), mu, sigma);
            return () -> Math.abs(s.sample());
        }
    }

    /**
     * Specialisation for the half-normal distribution.
     *
     * <p>Elimination of the {@code mu} location parameter simplifies the probability
     * functions and allows computation of the log density and inverse CDF/SF.
     */
    private static class HalfNormalDistribution extends FoldedNormalDistribution {
        /** Variance constant (1 - 2/pi). Computed using Matlab's VPA to 30 digits. */
        private static final double VAR = 0.36338022763241865692446494650994;
        /** The value of {@code log(sigma) + 0.5 * log(2*PI)} stored for faster computation. */
        private final double logSigmaPlusHalfLog2Pi;

        /**
         * @param sigma Scale parameter.
         */
        HalfNormalDistribution(double sigma) {
            super(sigma);
            logSigmaPlusHalfLog2Pi = Math.log(sigma) + Constants.HALF_LOG_TWO_PI;
        }

        @Override
        public double getMu() {
            return 0;
        }

        @Override
        public double density(double x) {
            if (x < 0) {
                return 0;
            }
            return 2 * ExtendedPrecision.expmhxx(x / sigma) / sigmaSqrt2pi;
        }

        @Override
        public double probability(double x0,
                                  double x1) {
            if (x0 > x1) {
                throw new DistributionException(DistributionException.INVALID_RANGE_LOW_GT_HIGH,
                                                x0, x1);
            }
            if (x0 <= 0) {
                return cumulativeProbability(x1);
            }
            // Assumes x1 >= x0 && x0 > 0
            return ErfDifference.value(x0 / sigmaSqrt2, x1 / sigmaSqrt2);
        }

        @Override
        public double logDensity(double x) {
            if (x < 0) {
                return Double.NEGATIVE_INFINITY;
            }
            final double z = x / sigma;
            return Constants.LN_TWO - 0.5 * z * z - logSigmaPlusHalfLog2Pi;
        }

        @Override
        public double cumulativeProbability(double x) {
            if (x <= 0) {
                return 0;
            }
            return Erf.value(x / sigmaSqrt2);
        }

        @Override
        public double survivalProbability(double x) {
            if (x <= 0) {
                return 1;
            }
            return Erfc.value(x / sigmaSqrt2);
        }

        @Override
        public double inverseCumulativeProbability(double p) {
            ArgumentUtils.checkProbability(p);
            // Addition of 0.0 ensures 0.0 is returned for p=-0.0
            return 0.0 + sigmaSqrt2 * InverseErf.value(p);
        }

        /** {@inheritDoc} */
        @Override
        public double inverseSurvivalProbability(double p) {
            ArgumentUtils.checkProbability(p);
            return sigmaSqrt2 * InverseErfc.value(p);
        }

        @Override
        public double getMean() {
            return sigma * Constants.ROOT_TWO_DIV_PI;
        }

        @Override
        public double getVariance() {
            // sigma^2 - mean^2
            // sigma^2 - (sigma^2 * 2/pi)
            return sigma * sigma * VAR;
        }

        @Override
        public Sampler createSampler(UniformRandomProvider rng) {
            // Return the absolute of a Gaussian distribution sampler.
            final SharedStateContinuousSampler s = ZigguratSampler.NormalizedGaussian.of(rng);
            return () -> Math.abs(s.sample() * sigma);
        }
    }

    /**
     * @param sigma Scale parameter.
     */
    FoldedNormalDistribution(double sigma) {
        this.sigma = sigma;
        // Minimise rounding error by computing sqrt(2 * sigma * sigma) exactly.
        // Compute using extended precision with care to avoid over/underflow.
        sigmaSqrt2 = ExtendedPrecision.sqrt2xx(sigma);
        // Compute sigma * sqrt(2 * pi)
        sigmaSqrt2pi = ExtendedPrecision.xsqrt2pi(sigma);
    }

    /**
     * Creates a folded normal distribution. If the location {@code mu} is zero this is
     * the half-normal distribution.
     *
     * @param mu Location parameter.
     * @param sigma Scale parameter.
     * @return the distribution
     * @throws IllegalArgumentException if {@code sigma <= 0}.
     */
    public static FoldedNormalDistribution of(double mu,
                                              double sigma) {
        if (sigma > 0) {
            if (mu == 0) {
                return new HalfNormalDistribution(sigma);
            }
            return new RegularFoldedNormalDistribution(mu, sigma);
        }
        // scale is zero, negative or nan
        throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, sigma);
    }

    /**
     * Gets the location parameter \( \mu \) of this distribution.
     *
     * @return the mu parameter.
     */
    public abstract double getMu();

    /**
     * Gets the scale parameter \( \sigma \) of this distribution.
     *
     * @return the sigma parameter.
     */
    public double getSigma() {
        return sigma;
    }

    /**
     * {@inheritDoc}
     *
     *
     * <p>For location parameter \( \mu \) and scale parameter \( \sigma \), the mean is:
     *
     * <p>\[ \sigma \sqrt{ \frac 2 \pi } \exp \left( \frac{-\mu^2}{2\sigma^2} \right) +
     *       \mu \operatorname{erf} \left( \frac \mu {\sqrt{2\sigma^2}} \right) \]
     *
     * <p>where \( \operatorname{erf} \) is the error function.
     */
    @Override
    public abstract double getMean();

    /**
     * {@inheritDoc}
     *
     * <p>For location parameter \( \mu \), scale parameter \( \sigma \) and a distribution
     * mean \( \mu_Y \), the variance is:
     *
     * <p>\[ \mu^2 + \sigma^2 - \mu_{Y}^2 \]
     */
    @Override
    public abstract double getVariance();

    /**
     * {@inheritDoc}
     *
     * <p>The lower bound of the support is always 0.
     *
     * @return 0.
     */
    @Override
    public double getSupportLowerBound() {
        return 0.0;
    }

    /**
     * {@inheritDoc}
     *
     * <p>The upper bound of the support is always positive infinity.
     *
     * @return {@linkplain Double#POSITIVE_INFINITY positive infinity}.
     */
    @Override
    public double getSupportUpperBound() {
        return Double.POSITIVE_INFINITY;
    }
}