You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2022/12/01 16:47:13 UTC
svn commit: r58392 [47/50] - in /dev/commons/statistics/1.0-RC1: ./ binaries/ site/ site/css/ site/images/ site/images/logos/ site/img/ site/js/ site/release-notes/ site/style/ site/userguide/ site/xref-test/ site/xref-test/org/ site/xref-test/org/apac...
Added: dev/commons/statistics/1.0-RC1/site/xref/org/apache/commons/statistics/distribution/TruncatedNormalDistribution.html
==============================================================================
--- dev/commons/statistics/1.0-RC1/site/xref/org/apache/commons/statistics/distribution/TruncatedNormalDistribution.html (added)
+++ dev/commons/statistics/1.0-RC1/site/xref/org/apache/commons/statistics/distribution/TruncatedNormalDistribution.html Thu Dec 1 16:47:12 2022
@@ -0,0 +1,596 @@
+<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
+<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
+<head><meta http-equiv="content-type" content="text/html; charset=UTF-8" />
+<title>TruncatedNormalDistribution xref</title>
+<link type="text/css" rel="stylesheet" href="../../../../../stylesheet.css" />
+</head>
+<body>
+<div id="overview"><a href="../../../../../../apidocs/org/apache/commons/statistics/distribution/TruncatedNormalDistribution.html">View Javadoc</a></div><pre>
+<a class="jxr_linenumber" name="L1" href="#L1">1</a> <em class="jxr_comment">/*</em>
+<a class="jxr_linenumber" name="L2" href="#L2">2</a> <em class="jxr_comment"> * Licensed to the Apache Software Foundation (ASF) under one or more</em>
+<a class="jxr_linenumber" name="L3" href="#L3">3</a> <em class="jxr_comment"> * contributor license agreements. See the NOTICE file distributed with</em>
+<a class="jxr_linenumber" name="L4" href="#L4">4</a> <em class="jxr_comment"> * this work for additional information regarding copyright ownership.</em>
+<a class="jxr_linenumber" name="L5" href="#L5">5</a> <em class="jxr_comment"> * The ASF licenses this file to You under the Apache License, Version 2.0</em>
+<a class="jxr_linenumber" name="L6" href="#L6">6</a> <em class="jxr_comment"> * (the "License"); you may not use this file except in compliance with</em>
+<a class="jxr_linenumber" name="L7" href="#L7">7</a> <em class="jxr_comment"> * the License. You may obtain a copy of the License at</em>
+<a class="jxr_linenumber" name="L8" href="#L8">8</a> <em class="jxr_comment"> *</em>
+<a class="jxr_linenumber" name="L9" href="#L9">9</a> <em class="jxr_comment"> * <a href="http://www.apache.org/licenses/LICENSE-2.0" target="alexandria_uri">http://www.apache.org/licenses/LICENSE-2.0</a></em>
+<a class="jxr_linenumber" name="L10" href="#L10">10</a> <em class="jxr_comment"> *</em>
+<a class="jxr_linenumber" name="L11" href="#L11">11</a> <em class="jxr_comment"> * Unless required by applicable law or agreed to in writing, software</em>
+<a class="jxr_linenumber" name="L12" href="#L12">12</a> <em class="jxr_comment"> * distributed under the License is distributed on an "AS IS" BASIS,</em>
+<a class="jxr_linenumber" name="L13" href="#L13">13</a> <em class="jxr_comment"> * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.</em>
+<a class="jxr_linenumber" name="L14" href="#L14">14</a> <em class="jxr_comment"> * See the License for the specific language governing permissions and</em>
+<a class="jxr_linenumber" name="L15" href="#L15">15</a> <em class="jxr_comment"> * limitations under the License.</em>
+<a class="jxr_linenumber" name="L16" href="#L16">16</a> <em class="jxr_comment"> */</em>
+<a class="jxr_linenumber" name="L17" href="#L17">17</a>
+<a class="jxr_linenumber" name="L18" href="#L18">18</a> <strong class="jxr_keyword">package</strong> org.apache.commons.statistics.distribution;
+<a class="jxr_linenumber" name="L19" href="#L19">19</a>
+<a class="jxr_linenumber" name="L20" href="#L20">20</a> <strong class="jxr_keyword">import</strong> java.util.function.DoubleSupplier;
+<a class="jxr_linenumber" name="L21" href="#L21">21</a> <strong class="jxr_keyword">import</strong> org.apache.commons.numbers.gamma.Erf;
+<a class="jxr_linenumber" name="L22" href="#L22">22</a> <strong class="jxr_keyword">import</strong> org.apache.commons.numbers.gamma.ErfDifference;
+<a class="jxr_linenumber" name="L23" href="#L23">23</a> <strong class="jxr_keyword">import</strong> org.apache.commons.numbers.gamma.Erfcx;
+<a class="jxr_linenumber" name="L24" href="#L24">24</a> <strong class="jxr_keyword">import</strong> org.apache.commons.rng.UniformRandomProvider;
+<a class="jxr_linenumber" name="L25" href="#L25">25</a> <strong class="jxr_keyword">import</strong> org.apache.commons.rng.sampling.distribution.ZigguratSampler;
+<a class="jxr_linenumber" name="L26" href="#L26">26</a>
+<a class="jxr_linenumber" name="L27" href="#L27">27</a> <em class="jxr_javadoccomment">/**</em>
+<a class="jxr_linenumber" name="L28" href="#L28">28</a> <em class="jxr_javadoccomment"> * Implementation of the truncated normal distribution.</em>
+<a class="jxr_linenumber" name="L29" href="#L29">29</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L30" href="#L30">30</a> <em class="jxr_javadoccomment"> * <p>The probability density function of \( X \) is:</em>
+<a class="jxr_linenumber" name="L31" href="#L31">31</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L32" href="#L32">32</a> <em class="jxr_javadoccomment"> * <p>\[ f(x;\mu,\sigma,a,b) = \frac{1}{\sigma}\,\frac{\phi(\frac{x - \mu}{\sigma})}{\Phi(\frac{b - \mu}{\sigma}) - \Phi(\frac{a - \mu}{\sigma}) } \]</em>
+<a class="jxr_linenumber" name="L33" href="#L33">33</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L34" href="#L34">34</a> <em class="jxr_javadoccomment"> * <p>for \( \mu \) mean of the parent normal distribution,</em>
+<a class="jxr_linenumber" name="L35" href="#L35">35</a> <em class="jxr_javadoccomment"> * \( \sigma \) standard deviation of the parent normal distribution,</em>
+<a class="jxr_linenumber" name="L36" href="#L36">36</a> <em class="jxr_javadoccomment"> * \( -\infty \le a \lt b \le \infty \) the truncation interval, and</em>
+<a class="jxr_linenumber" name="L37" href="#L37">37</a> <em class="jxr_javadoccomment"> * \( x \in [a, b] \), where \( \phi \) is the probability</em>
+<a class="jxr_linenumber" name="L38" href="#L38">38</a> <em class="jxr_javadoccomment"> * density function of the standard normal distribution and \( \Phi \)</em>
+<a class="jxr_linenumber" name="L39" href="#L39">39</a> <em class="jxr_javadoccomment"> * is its cumulative distribution function.</em>
+<a class="jxr_linenumber" name="L40" href="#L40">40</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L41" href="#L41">41</a> <em class="jxr_javadoccomment"> * @see <a href="<a href="https://en.wikipedia.org/wiki/Truncated_normal_distribution" target="alexandria_uri">https://en.wikipedia.org/wiki/Truncated_normal_distribution</a>"></em>
+<a class="jxr_linenumber" name="L42" href="#L42">42</a> <em class="jxr_javadoccomment"> * Truncated normal distribution (Wikipedia)</a></em>
+<a class="jxr_linenumber" name="L43" href="#L43">43</a> <em class="jxr_javadoccomment"> */</em>
+<a class="jxr_linenumber" name="L44" href="#L44">44</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">class</strong> <a name="TruncatedNormalDistribution" href="../../../../../org/apache/commons/statistics/distribution/TruncatedNormalDistribution.html#TruncatedNormalDistribution">TruncatedNormalDistribution</a> <strong class="jxr_keyword">extends</strong> <a name="AbstractContinuousDistribution" href="../../../../../org/apache/commons/statistics/distribution/AbstractContinuousDistribution.html#AbstractContinuousDistribution">AbstractContinuousDistribution</a> {
+<a class="jxr_linenumber" name="L45" href="#L45">45</a>
+<a class="jxr_linenumber" name="L46" href="#L46">46</a> <em class="jxr_javadoccomment">/** The max allowed value for x where (x*x) will not overflow.</em>
+<a class="jxr_linenumber" name="L47" href="#L47">47</a> <em class="jxr_javadoccomment"> * This is a limit on computation of the moments of the truncated normal</em>
+<a class="jxr_linenumber" name="L48" href="#L48">48</a> <em class="jxr_javadoccomment"> * as some calculations assume x*x is finite. Value is sqrt(MAX_VALUE). */</em>
+<a class="jxr_linenumber" name="L49" href="#L49">49</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> MAX_X = 0x1.fffffffffffffp511;
+<a class="jxr_linenumber" name="L50" href="#L50">50</a>
+<a class="jxr_linenumber" name="L51" href="#L51">51</a> <em class="jxr_javadoccomment">/** The min allowed probability range of the parent normal distribution.</em>
+<a class="jxr_linenumber" name="L52" href="#L52">52</a> <em class="jxr_javadoccomment"> * Set to 0.0. This may be too low for accurate usage. It is a signal that</em>
+<a class="jxr_linenumber" name="L53" href="#L53">53</a> <em class="jxr_javadoccomment"> * the truncation is invalid. */</em>
+<a class="jxr_linenumber" name="L54" href="#L54">54</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> MIN_P = 0.0;
+<a class="jxr_linenumber" name="L55" href="#L55">55</a>
+<a class="jxr_linenumber" name="L56" href="#L56">56</a> <em class="jxr_javadoccomment">/** sqrt(2). */</em>
+<a class="jxr_linenumber" name="L57" href="#L57">57</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> ROOT2 = 1.414213562373095048801688724209698078;
+<a class="jxr_linenumber" name="L58" href="#L58">58</a> <em class="jxr_javadoccomment">/** Normalisation constant 2 / sqrt(2 pi) = sqrt(2 / pi). */</em>
+<a class="jxr_linenumber" name="L59" href="#L59">59</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> ROOT_2_PI = 0.797884560802865405726436165423365309;
+<a class="jxr_linenumber" name="L60" href="#L60">60</a> <em class="jxr_javadoccomment">/** Normalisation constant sqrt(2 pi) / 2 = sqrt(pi / 2). */</em>
+<a class="jxr_linenumber" name="L61" href="#L61">61</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> ROOT_PI_2 = 1.253314137315500251207882642405522626;
+<a class="jxr_linenumber" name="L62" href="#L62">62</a>
+<a class="jxr_linenumber" name="L63" href="#L63">63</a> <em class="jxr_javadoccomment">/**</em>
+<a class="jxr_linenumber" name="L64" href="#L64">64</a> <em class="jxr_javadoccomment"> * The threshold to switch to a rejection sampler. When the truncated</em>
+<a class="jxr_linenumber" name="L65" href="#L65">65</a> <em class="jxr_javadoccomment"> * distribution covers more than this fraction of the CDF then rejection</em>
+<a class="jxr_linenumber" name="L66" href="#L66">66</a> <em class="jxr_javadoccomment"> * sampling will be more efficient than inverse CDF sampling. Performance</em>
+<a class="jxr_linenumber" name="L67" href="#L67">67</a> <em class="jxr_javadoccomment"> * benchmarks indicate that a normalized Gaussian sampler is up to 10 times</em>
+<a class="jxr_linenumber" name="L68" href="#L68">68</a> <em class="jxr_javadoccomment"> * faster than inverse transform sampling using a fast random generator. See</em>
+<a class="jxr_linenumber" name="L69" href="#L69">69</a> <em class="jxr_javadoccomment"> * STATISTICS-55.</em>
+<a class="jxr_linenumber" name="L70" href="#L70">70</a> <em class="jxr_javadoccomment"> */</em>
+<a class="jxr_linenumber" name="L71" href="#L71">71</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> REJECTION_THRESHOLD = 0.2;
+<a class="jxr_linenumber" name="L72" href="#L72">72</a>
+<a class="jxr_linenumber" name="L73" href="#L73">73</a> <em class="jxr_javadoccomment">/** Parent normal distribution. */</em>
+<a class="jxr_linenumber" name="L74" href="#L74">74</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">final</strong> <a name="NormalDistribution" href="../../../../../org/apache/commons/statistics/distribution/NormalDistribution.html#NormalDistribution">NormalDistribution</a> parentNormal;
+<a class="jxr_linenumber" name="L75" href="#L75">75</a> <em class="jxr_javadoccomment">/** Lower bound of this distribution. */</em>
+<a class="jxr_linenumber" name="L76" href="#L76">76</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> lower;
+<a class="jxr_linenumber" name="L77" href="#L77">77</a> <em class="jxr_javadoccomment">/** Upper bound of this distribution. */</em>
+<a class="jxr_linenumber" name="L78" href="#L78">78</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> upper;
+<a class="jxr_linenumber" name="L79" href="#L79">79</a>
+<a class="jxr_linenumber" name="L80" href="#L80">80</a> <em class="jxr_javadoccomment">/** Stored value of {@code parentNormal.probability(lower, upper)}. This is used to</em>
+<a class="jxr_linenumber" name="L81" href="#L81">81</a> <em class="jxr_javadoccomment"> * normalise the probability computations. */</em>
+<a class="jxr_linenumber" name="L82" href="#L82">82</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> cdfDelta;
+<a class="jxr_linenumber" name="L83" href="#L83">83</a> <em class="jxr_javadoccomment">/** log(cdfDelta). */</em>
+<a class="jxr_linenumber" name="L84" href="#L84">84</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> logCdfDelta;
+<a class="jxr_linenumber" name="L85" href="#L85">85</a> <em class="jxr_javadoccomment">/** Stored value of {@code parentNormal.cumulativeProbability(lower)}. Used to map</em>
+<a class="jxr_linenumber" name="L86" href="#L86">86</a> <em class="jxr_javadoccomment"> * a probability into the range of the parent normal distribution. */</em>
+<a class="jxr_linenumber" name="L87" href="#L87">87</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> cdfAlpha;
+<a class="jxr_linenumber" name="L88" href="#L88">88</a> <em class="jxr_javadoccomment">/** Stored value of {@code parentNormal.survivalProbability(upper)}. Used to map</em>
+<a class="jxr_linenumber" name="L89" href="#L89">89</a> <em class="jxr_javadoccomment"> * a probability into the range of the parent normal distribution. */</em>
+<a class="jxr_linenumber" name="L90" href="#L90">90</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> sfBeta;
+<a class="jxr_linenumber" name="L91" href="#L91">91</a>
+<a class="jxr_linenumber" name="L92" href="#L92">92</a> <em class="jxr_javadoccomment">/**</em>
+<a class="jxr_linenumber" name="L93" href="#L93">93</a> <em class="jxr_javadoccomment"> * @param parent Parent distribution.</em>
+<a class="jxr_linenumber" name="L94" href="#L94">94</a> <em class="jxr_javadoccomment"> * @param z Probability of the parent distribution for {@code [lower, upper]}.</em>
+<a class="jxr_linenumber" name="L95" href="#L95">95</a> <em class="jxr_javadoccomment"> * @param lower Lower bound (inclusive) of the distribution, can be {@link Double#NEGATIVE_INFINITY}.</em>
+<a class="jxr_linenumber" name="L96" href="#L96">96</a> <em class="jxr_javadoccomment"> * @param upper Upper bound (inclusive) of the distribution, can be {@link Double#POSITIVE_INFINITY}.</em>
+<a class="jxr_linenumber" name="L97" href="#L97">97</a> <em class="jxr_javadoccomment"> */</em>
+<a class="jxr_linenumber" name="L98" href="#L98">98</a> <strong class="jxr_keyword">private</strong> <a name="TruncatedNormalDistribution" href="../../../../../org/apache/commons/statistics/distribution/TruncatedNormalDistribution.html#TruncatedNormalDistribution">TruncatedNormalDistribution</a>(<a name="NormalDistribution" href="../../../../../org/apache/commons/statistics/distribution/NormalDistribution.html#NormalDistribution">NormalDistribution</a> parent, <strong class="jxr_keyword">double</strong> z, <strong class="jxr_keyword">double</strong> lower, <strong class="jxr_keyword">double</strong> upper) {
+<a class="jxr_linenumber" name="L99" href="#L99">99</a> <strong class="jxr_keyword">this</strong>.parentNormal = parent;
+<a class="jxr_linenumber" name="L100" href="#L100">100</a> <strong class="jxr_keyword">this</strong>.lower = lower;
+<a class="jxr_linenumber" name="L101" href="#L101">101</a> <strong class="jxr_keyword">this</strong>.upper = upper;
+<a class="jxr_linenumber" name="L102" href="#L102">102</a>
+<a class="jxr_linenumber" name="L103" href="#L103">103</a> cdfDelta = z;
+<a class="jxr_linenumber" name="L104" href="#L104">104</a> logCdfDelta = Math.log(cdfDelta);
+<a class="jxr_linenumber" name="L105" href="#L105">105</a> <em class="jxr_comment">// Used to map the inverse probability.</em>
+<a class="jxr_linenumber" name="L106" href="#L106">106</a> cdfAlpha = parentNormal.cumulativeProbability(lower);
+<a class="jxr_linenumber" name="L107" href="#L107">107</a> sfBeta = parentNormal.survivalProbability(upper);
+<a class="jxr_linenumber" name="L108" href="#L108">108</a> }
+<a class="jxr_linenumber" name="L109" href="#L109">109</a>
+<a class="jxr_linenumber" name="L110" href="#L110">110</a> <em class="jxr_javadoccomment">/**</em>
+<a class="jxr_linenumber" name="L111" href="#L111">111</a> <em class="jxr_javadoccomment"> * Creates a truncated normal distribution.</em>
+<a class="jxr_linenumber" name="L112" href="#L112">112</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L113" href="#L113">113</a> <em class="jxr_javadoccomment"> * <p>Note that the {@code mean} and {@code sd} is of the parent normal distribution,</em>
+<a class="jxr_linenumber" name="L114" href="#L114">114</a> <em class="jxr_javadoccomment"> * and not the true mean and standard deviation of the truncated normal distribution.</em>
+<a class="jxr_linenumber" name="L115" href="#L115">115</a> <em class="jxr_javadoccomment"> * The {@code lower} and {@code upper} bounds define the truncation of the parent</em>
+<a class="jxr_linenumber" name="L116" href="#L116">116</a> <em class="jxr_javadoccomment"> * normal distribution.</em>
+<a class="jxr_linenumber" name="L117" href="#L117">117</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L118" href="#L118">118</a> <em class="jxr_javadoccomment"> * @param mean Mean for the parent distribution.</em>
+<a class="jxr_linenumber" name="L119" href="#L119">119</a> <em class="jxr_javadoccomment"> * @param sd Standard deviation for the parent distribution.</em>
+<a class="jxr_linenumber" name="L120" href="#L120">120</a> <em class="jxr_javadoccomment"> * @param lower Lower bound (inclusive) of the distribution, can be {@link Double#NEGATIVE_INFINITY}.</em>
+<a class="jxr_linenumber" name="L121" href="#L121">121</a> <em class="jxr_javadoccomment"> * @param upper Upper bound (inclusive) of the distribution, can be {@link Double#POSITIVE_INFINITY}.</em>
+<a class="jxr_linenumber" name="L122" href="#L122">122</a> <em class="jxr_javadoccomment"> * @return the distribution</em>
+<a class="jxr_linenumber" name="L123" href="#L123">123</a> <em class="jxr_javadoccomment"> * @throws IllegalArgumentException if {@code sd <= 0}; if {@code lower >= upper}; or if</em>
+<a class="jxr_linenumber" name="L124" href="#L124">124</a> <em class="jxr_javadoccomment"> * the truncation covers no probability range in the parent distribution.</em>
+<a class="jxr_linenumber" name="L125" href="#L125">125</a> <em class="jxr_javadoccomment"> */</em>
+<a class="jxr_linenumber" name="L126" href="#L126">126</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">static</strong> <a name="TruncatedNormalDistribution" href="../../../../../org/apache/commons/statistics/distribution/TruncatedNormalDistribution.html#TruncatedNormalDistribution">TruncatedNormalDistribution</a> of(<strong class="jxr_keyword">double</strong> mean, <strong class="jxr_keyword">double</strong> sd, <strong class="jxr_keyword">double</strong> lower, <strong class="jxr_keyword">double</strong> upper) {
+<a class="jxr_linenumber" name="L127" href="#L127">127</a> <strong class="jxr_keyword">if</strong> (sd <= 0) {
+<a class="jxr_linenumber" name="L128" href="#L128">128</a> <strong class="jxr_keyword">throw</strong> <strong class="jxr_keyword">new</strong> <a name="DistributionException" href="../../../../../org/apache/commons/statistics/distribution/DistributionException.html#DistributionException">DistributionException</a>(DistributionException.NOT_STRICTLY_POSITIVE, sd);
+<a class="jxr_linenumber" name="L129" href="#L129">129</a> }
+<a class="jxr_linenumber" name="L130" href="#L130">130</a> <strong class="jxr_keyword">if</strong> (lower >= upper) {
+<a class="jxr_linenumber" name="L131" href="#L131">131</a> <strong class="jxr_keyword">throw</strong> <strong class="jxr_keyword">new</strong> <a name="DistributionException" href="../../../../../org/apache/commons/statistics/distribution/DistributionException.html#DistributionException">DistributionException</a>(DistributionException.INVALID_RANGE_LOW_GTE_HIGH, lower, upper);
+<a class="jxr_linenumber" name="L132" href="#L132">132</a> }
+<a class="jxr_linenumber" name="L133" href="#L133">133</a>
+<a class="jxr_linenumber" name="L134" href="#L134">134</a> <em class="jxr_comment">// Use an instance for the parent normal distribution to maximise accuracy</em>
+<a class="jxr_linenumber" name="L135" href="#L135">135</a> <em class="jxr_comment">// in range computations using the error function</em>
+<a class="jxr_linenumber" name="L136" href="#L136">136</a> <strong class="jxr_keyword">final</strong> <a name="NormalDistribution" href="../../../../../org/apache/commons/statistics/distribution/NormalDistribution.html#NormalDistribution">NormalDistribution</a> parent = NormalDistribution.of(mean, sd);
+<a class="jxr_linenumber" name="L137" href="#L137">137</a>
+<a class="jxr_linenumber" name="L138" href="#L138">138</a> <em class="jxr_comment">// If there is no computable range then raise an exception.</em>
+<a class="jxr_linenumber" name="L139" href="#L139">139</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> z = parent.probability(lower, upper);
+<a class="jxr_linenumber" name="L140" href="#L140">140</a> <strong class="jxr_keyword">if</strong> (z <= MIN_P) {
+<a class="jxr_linenumber" name="L141" href="#L141">141</a> <em class="jxr_comment">// Map the bounds to a standard normal distribution for the message</em>
+<a class="jxr_linenumber" name="L142" href="#L142">142</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> a = (lower - mean) / sd;
+<a class="jxr_linenumber" name="L143" href="#L143">143</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> b = (upper - mean) / sd;
+<a class="jxr_linenumber" name="L144" href="#L144">144</a> <strong class="jxr_keyword">throw</strong> <strong class="jxr_keyword">new</strong> <a name="DistributionException" href="../../../../../org/apache/commons/statistics/distribution/DistributionException.html#DistributionException">DistributionException</a>(
+<a class="jxr_linenumber" name="L145" href="#L145">145</a> <span class="jxr_string">"Excess truncation of standard normal : CDF(%s, %s) = %s"</span>, a, b, z);
+<a class="jxr_linenumber" name="L146" href="#L146">146</a> }
+<a class="jxr_linenumber" name="L147" href="#L147">147</a>
+<a class="jxr_linenumber" name="L148" href="#L148">148</a> <em class="jxr_comment">// Here we have a meaningful truncation. Note that excess truncation may not be optimal.</em>
+<a class="jxr_linenumber" name="L149" href="#L149">149</a> <em class="jxr_comment">// For example truncation close to zero where the PDF is constant can be approximated</em>
+<a class="jxr_linenumber" name="L150" href="#L150">150</a> <em class="jxr_comment">// using a uniform distribution.</em>
+<a class="jxr_linenumber" name="L151" href="#L151">151</a>
+<a class="jxr_linenumber" name="L152" href="#L152">152</a> <strong class="jxr_keyword">return</strong> <strong class="jxr_keyword">new</strong> <a name="TruncatedNormalDistribution" href="../../../../../org/apache/commons/statistics/distribution/TruncatedNormalDistribution.html#TruncatedNormalDistribution">TruncatedNormalDistribution</a>(parent, z, lower, upper);
+<a class="jxr_linenumber" name="L153" href="#L153">153</a> }
+<a class="jxr_linenumber" name="L154" href="#L154">154</a>
+<a class="jxr_linenumber" name="L155" href="#L155">155</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em>
+<a class="jxr_linenumber" name="L156" href="#L156">156</a> @Override
+<a class="jxr_linenumber" name="L157" href="#L157">157</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> density(<strong class="jxr_keyword">double</strong> x) {
+<a class="jxr_linenumber" name="L158" href="#L158">158</a> <strong class="jxr_keyword">if</strong> (x < lower || x > upper) {
+<a class="jxr_linenumber" name="L159" href="#L159">159</a> <strong class="jxr_keyword">return</strong> 0;
+<a class="jxr_linenumber" name="L160" href="#L160">160</a> }
+<a class="jxr_linenumber" name="L161" href="#L161">161</a> <strong class="jxr_keyword">return</strong> parentNormal.density(x) / cdfDelta;
+<a class="jxr_linenumber" name="L162" href="#L162">162</a> }
+<a class="jxr_linenumber" name="L163" href="#L163">163</a>
+<a class="jxr_linenumber" name="L164" href="#L164">164</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em>
+<a class="jxr_linenumber" name="L165" href="#L165">165</a> @Override
+<a class="jxr_linenumber" name="L166" href="#L166">166</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> probability(<strong class="jxr_keyword">double</strong> x0, <strong class="jxr_keyword">double</strong> x1) {
+<a class="jxr_linenumber" name="L167" href="#L167">167</a> <strong class="jxr_keyword">if</strong> (x0 > x1) {
+<a class="jxr_linenumber" name="L168" href="#L168">168</a> <strong class="jxr_keyword">throw</strong> <strong class="jxr_keyword">new</strong> <a name="DistributionException" href="../../../../../org/apache/commons/statistics/distribution/DistributionException.html#DistributionException">DistributionException</a>(DistributionException.INVALID_RANGE_LOW_GT_HIGH,
+<a class="jxr_linenumber" name="L169" href="#L169">169</a> x0, x1);
+<a class="jxr_linenumber" name="L170" href="#L170">170</a> }
+<a class="jxr_linenumber" name="L171" href="#L171">171</a> <strong class="jxr_keyword">return</strong> parentNormal.probability(clipToRange(x0), clipToRange(x1)) / cdfDelta;
+<a class="jxr_linenumber" name="L172" href="#L172">172</a> }
+<a class="jxr_linenumber" name="L173" href="#L173">173</a>
+<a class="jxr_linenumber" name="L174" href="#L174">174</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em>
+<a class="jxr_linenumber" name="L175" href="#L175">175</a> @Override
+<a class="jxr_linenumber" name="L176" href="#L176">176</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> logDensity(<strong class="jxr_keyword">double</strong> x) {
+<a class="jxr_linenumber" name="L177" href="#L177">177</a> <strong class="jxr_keyword">if</strong> (x < lower || x > upper) {
+<a class="jxr_linenumber" name="L178" href="#L178">178</a> <strong class="jxr_keyword">return</strong> Double.NEGATIVE_INFINITY;
+<a class="jxr_linenumber" name="L179" href="#L179">179</a> }
+<a class="jxr_linenumber" name="L180" href="#L180">180</a> <strong class="jxr_keyword">return</strong> parentNormal.logDensity(x) - logCdfDelta;
+<a class="jxr_linenumber" name="L181" href="#L181">181</a> }
+<a class="jxr_linenumber" name="L182" href="#L182">182</a>
+<a class="jxr_linenumber" name="L183" href="#L183">183</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em>
+<a class="jxr_linenumber" name="L184" href="#L184">184</a> @Override
+<a class="jxr_linenumber" name="L185" href="#L185">185</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> cumulativeProbability(<strong class="jxr_keyword">double</strong> x) {
+<a class="jxr_linenumber" name="L186" href="#L186">186</a> <strong class="jxr_keyword">if</strong> (x <= lower) {
+<a class="jxr_linenumber" name="L187" href="#L187">187</a> <strong class="jxr_keyword">return</strong> 0;
+<a class="jxr_linenumber" name="L188" href="#L188">188</a> } <strong class="jxr_keyword">else</strong> <strong class="jxr_keyword">if</strong> (x >= upper) {
+<a class="jxr_linenumber" name="L189" href="#L189">189</a> <strong class="jxr_keyword">return</strong> 1;
+<a class="jxr_linenumber" name="L190" href="#L190">190</a> }
+<a class="jxr_linenumber" name="L191" href="#L191">191</a> <strong class="jxr_keyword">return</strong> parentNormal.probability(lower, x) / cdfDelta;
+<a class="jxr_linenumber" name="L192" href="#L192">192</a> }
+<a class="jxr_linenumber" name="L193" href="#L193">193</a>
+<a class="jxr_linenumber" name="L194" href="#L194">194</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em>
+<a class="jxr_linenumber" name="L195" href="#L195">195</a> @Override
+<a class="jxr_linenumber" name="L196" href="#L196">196</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> survivalProbability(<strong class="jxr_keyword">double</strong> x) {
+<a class="jxr_linenumber" name="L197" href="#L197">197</a> <strong class="jxr_keyword">if</strong> (x <= lower) {
+<a class="jxr_linenumber" name="L198" href="#L198">198</a> <strong class="jxr_keyword">return</strong> 1;
+<a class="jxr_linenumber" name="L199" href="#L199">199</a> } <strong class="jxr_keyword">else</strong> <strong class="jxr_keyword">if</strong> (x >= upper) {
+<a class="jxr_linenumber" name="L200" href="#L200">200</a> <strong class="jxr_keyword">return</strong> 0;
+<a class="jxr_linenumber" name="L201" href="#L201">201</a> }
+<a class="jxr_linenumber" name="L202" href="#L202">202</a> <strong class="jxr_keyword">return</strong> parentNormal.probability(x, upper) / cdfDelta;
+<a class="jxr_linenumber" name="L203" href="#L203">203</a> }
+<a class="jxr_linenumber" name="L204" href="#L204">204</a>
+<a class="jxr_linenumber" name="L205" href="#L205">205</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em>
+<a class="jxr_linenumber" name="L206" href="#L206">206</a> @Override
+<a class="jxr_linenumber" name="L207" href="#L207">207</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> inverseCumulativeProbability(<strong class="jxr_keyword">double</strong> p) {
+<a class="jxr_linenumber" name="L208" href="#L208">208</a> ArgumentUtils.checkProbability(p);
+<a class="jxr_linenumber" name="L209" href="#L209">209</a> <em class="jxr_comment">// Exact bound</em>
+<a class="jxr_linenumber" name="L210" href="#L210">210</a> <strong class="jxr_keyword">if</strong> (p == 0) {
+<a class="jxr_linenumber" name="L211" href="#L211">211</a> <strong class="jxr_keyword">return</strong> lower;
+<a class="jxr_linenumber" name="L212" href="#L212">212</a> } <strong class="jxr_keyword">else</strong> <strong class="jxr_keyword">if</strong> (p == 1) {
+<a class="jxr_linenumber" name="L213" href="#L213">213</a> <strong class="jxr_keyword">return</strong> upper;
+<a class="jxr_linenumber" name="L214" href="#L214">214</a> }
+<a class="jxr_linenumber" name="L215" href="#L215">215</a> <em class="jxr_comment">// Linearly map p to the range [lower, upper]</em>
+<a class="jxr_linenumber" name="L216" href="#L216">216</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> x = parentNormal.inverseCumulativeProbability(cdfAlpha + p * cdfDelta);
+<a class="jxr_linenumber" name="L217" href="#L217">217</a> <strong class="jxr_keyword">return</strong> clipToRange(x);
+<a class="jxr_linenumber" name="L218" href="#L218">218</a> }
+<a class="jxr_linenumber" name="L219" href="#L219">219</a>
+<a class="jxr_linenumber" name="L220" href="#L220">220</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em>
+<a class="jxr_linenumber" name="L221" href="#L221">221</a> @Override
+<a class="jxr_linenumber" name="L222" href="#L222">222</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> inverseSurvivalProbability(<strong class="jxr_keyword">double</strong> p) {
+<a class="jxr_linenumber" name="L223" href="#L223">223</a> ArgumentUtils.checkProbability(p);
+<a class="jxr_linenumber" name="L224" href="#L224">224</a> <em class="jxr_comment">// Exact bound</em>
+<a class="jxr_linenumber" name="L225" href="#L225">225</a> <strong class="jxr_keyword">if</strong> (p == 1) {
+<a class="jxr_linenumber" name="L226" href="#L226">226</a> <strong class="jxr_keyword">return</strong> lower;
+<a class="jxr_linenumber" name="L227" href="#L227">227</a> } <strong class="jxr_keyword">else</strong> <strong class="jxr_keyword">if</strong> (p == 0) {
+<a class="jxr_linenumber" name="L228" href="#L228">228</a> <strong class="jxr_keyword">return</strong> upper;
+<a class="jxr_linenumber" name="L229" href="#L229">229</a> }
+<a class="jxr_linenumber" name="L230" href="#L230">230</a> <em class="jxr_comment">// Linearly map p to the range [lower, upper]</em>
+<a class="jxr_linenumber" name="L231" href="#L231">231</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> x = parentNormal.inverseSurvivalProbability(sfBeta + p * cdfDelta);
+<a class="jxr_linenumber" name="L232" href="#L232">232</a> <strong class="jxr_keyword">return</strong> clipToRange(x);
+<a class="jxr_linenumber" name="L233" href="#L233">233</a> }
+<a class="jxr_linenumber" name="L234" href="#L234">234</a>
+<a class="jxr_linenumber" name="L235" href="#L235">235</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em>
+<a class="jxr_linenumber" name="L236" href="#L236">236</a> @Override
+<a class="jxr_linenumber" name="L237" href="#L237">237</a> <strong class="jxr_keyword">public</strong> Sampler createSampler(UniformRandomProvider rng) {
+<a class="jxr_linenumber" name="L238" href="#L238">238</a> <em class="jxr_comment">// If the truncation covers a reasonable amount of the normal distribution</em>
+<a class="jxr_linenumber" name="L239" href="#L239">239</a> <em class="jxr_comment">// then a rejection sampler can be used.</em>
+<a class="jxr_linenumber" name="L240" href="#L240">240</a> <strong class="jxr_keyword">double</strong> threshold = REJECTION_THRESHOLD;
+<a class="jxr_linenumber" name="L241" href="#L241">241</a> <em class="jxr_comment">// If the truncation is entirely in the upper or lower half then adjust the</em>
+<a class="jxr_linenumber" name="L242" href="#L242">242</a> <em class="jxr_comment">// threshold as twice the samples can be used</em>
+<a class="jxr_linenumber" name="L243" href="#L243">243</a> <strong class="jxr_keyword">if</strong> (lower >= 0 || upper <= 0) {
+<a class="jxr_linenumber" name="L244" href="#L244">244</a> threshold *= 0.5;
+<a class="jxr_linenumber" name="L245" href="#L245">245</a> }
+<a class="jxr_linenumber" name="L246" href="#L246">246</a>
+<a class="jxr_linenumber" name="L247" href="#L247">247</a> <strong class="jxr_keyword">if</strong> (cdfDelta > threshold) {
+<a class="jxr_linenumber" name="L248" href="#L248">248</a> <em class="jxr_comment">// Create the rejection sampler</em>
+<a class="jxr_linenumber" name="L249" href="#L249">249</a> <strong class="jxr_keyword">final</strong> ZigguratSampler.NormalizedGaussian sampler = ZigguratSampler.NormalizedGaussian.of(rng);
+<a class="jxr_linenumber" name="L250" href="#L250">250</a> DoubleSupplier gen;
+<a class="jxr_linenumber" name="L251" href="#L251">251</a> <em class="jxr_comment">// Use mirroring if possible</em>
+<a class="jxr_linenumber" name="L252" href="#L252">252</a> <strong class="jxr_keyword">if</strong> (lower >= 0) {
+<a class="jxr_linenumber" name="L253" href="#L253">253</a> <em class="jxr_comment">// Return the upper-half of the Gaussian</em>
+<a class="jxr_linenumber" name="L254" href="#L254">254</a> gen = () -> Math.abs(sampler.sample());
+<a class="jxr_linenumber" name="L255" href="#L255">255</a> } <strong class="jxr_keyword">else</strong> <strong class="jxr_keyword">if</strong> (upper <= 0) {
+<a class="jxr_linenumber" name="L256" href="#L256">256</a> <em class="jxr_comment">// Return the lower-half of the Gaussian</em>
+<a class="jxr_linenumber" name="L257" href="#L257">257</a> gen = () -> -Math.abs(sampler.sample());
+<a class="jxr_linenumber" name="L258" href="#L258">258</a> } <strong class="jxr_keyword">else</strong> {
+<a class="jxr_linenumber" name="L259" href="#L259">259</a> <em class="jxr_comment">// Return the full range of the Gaussian</em>
+<a class="jxr_linenumber" name="L260" href="#L260">260</a> gen = sampler::sample;
+<a class="jxr_linenumber" name="L261" href="#L261">261</a> }
+<a class="jxr_linenumber" name="L262" href="#L262">262</a> <em class="jxr_comment">// Map the bounds to a standard normal distribution</em>
+<a class="jxr_linenumber" name="L263" href="#L263">263</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> u = parentNormal.getMean();
+<a class="jxr_linenumber" name="L264" href="#L264">264</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> s = parentNormal.getStandardDeviation();
+<a class="jxr_linenumber" name="L265" href="#L265">265</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> a = (lower - u) / s;
+<a class="jxr_linenumber" name="L266" href="#L266">266</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> b = (upper - u) / s;
+<a class="jxr_linenumber" name="L267" href="#L267">267</a> <em class="jxr_comment">// Sample in [a, b] using rejection</em>
+<a class="jxr_linenumber" name="L268" href="#L268">268</a> <strong class="jxr_keyword">return</strong> () -> {
+<a class="jxr_linenumber" name="L269" href="#L269">269</a> <strong class="jxr_keyword">double</strong> x = gen.getAsDouble();
+<a class="jxr_linenumber" name="L270" href="#L270">270</a> <strong class="jxr_keyword">while</strong> (x < a || x > b) {
+<a class="jxr_linenumber" name="L271" href="#L271">271</a> x = gen.getAsDouble();
+<a class="jxr_linenumber" name="L272" href="#L272">272</a> }
+<a class="jxr_linenumber" name="L273" href="#L273">273</a> <em class="jxr_comment">// Avoid floating-point error when mapping back</em>
+<a class="jxr_linenumber" name="L274" href="#L274">274</a> <strong class="jxr_keyword">return</strong> clipToRange(u + x * s);
+<a class="jxr_linenumber" name="L275" href="#L275">275</a> };
+<a class="jxr_linenumber" name="L276" href="#L276">276</a> }
+<a class="jxr_linenumber" name="L277" href="#L277">277</a>
+<a class="jxr_linenumber" name="L278" href="#L278">278</a> <em class="jxr_comment">// Default to an inverse CDF sampler</em>
+<a class="jxr_linenumber" name="L279" href="#L279">279</a> <strong class="jxr_keyword">return</strong> <strong class="jxr_keyword">super</strong>.createSampler(rng);
+<a class="jxr_linenumber" name="L280" href="#L280">280</a> }
+<a class="jxr_linenumber" name="L281" href="#L281">281</a>
+<a class="jxr_linenumber" name="L282" href="#L282">282</a> <em class="jxr_javadoccomment">/**</em>
+<a class="jxr_linenumber" name="L283" href="#L283">283</a> <em class="jxr_javadoccomment"> * {@inheritDoc}</em>
+<a class="jxr_linenumber" name="L284" href="#L284">284</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L285" href="#L285">285</a> <em class="jxr_javadoccomment"> * <p>Represents the true mean of the truncated normal distribution rather</em>
+<a class="jxr_linenumber" name="L286" href="#L286">286</a> <em class="jxr_javadoccomment"> * than the parent normal distribution mean.</em>
+<a class="jxr_linenumber" name="L287" href="#L287">287</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L288" href="#L288">288</a> <em class="jxr_javadoccomment"> * <p>For \( \mu \) mean of the parent normal distribution,</em>
+<a class="jxr_linenumber" name="L289" href="#L289">289</a> <em class="jxr_javadoccomment"> * \( \sigma \) standard deviation of the parent normal distribution, and</em>
+<a class="jxr_linenumber" name="L290" href="#L290">290</a> <em class="jxr_javadoccomment"> * \( a \lt b \) the truncation interval of the parent normal distribution, the mean is:</em>
+<a class="jxr_linenumber" name="L291" href="#L291">291</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L292" href="#L292">292</a> <em class="jxr_javadoccomment"> * <p>\[ \mu + \frac{\phi(a)-\phi(b)}{\Phi(b) - \Phi(a)}\sigma \]</em>
+<a class="jxr_linenumber" name="L293" href="#L293">293</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L294" href="#L294">294</a> <em class="jxr_javadoccomment"> * <p>where \( \phi \) is the probability density function of the standard normal distribution</em>
+<a class="jxr_linenumber" name="L295" href="#L295">295</a> <em class="jxr_javadoccomment"> * and \( \Phi \) is its cumulative distribution function.</em>
+<a class="jxr_linenumber" name="L296" href="#L296">296</a> <em class="jxr_javadoccomment"> */</em>
+<a class="jxr_linenumber" name="L297" href="#L297">297</a> @Override
+<a class="jxr_linenumber" name="L298" href="#L298">298</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> getMean() {
+<a class="jxr_linenumber" name="L299" href="#L299">299</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> u = parentNormal.getMean();
+<a class="jxr_linenumber" name="L300" href="#L300">300</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> s = parentNormal.getStandardDeviation();
+<a class="jxr_linenumber" name="L301" href="#L301">301</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> a = (lower - u) / s;
+<a class="jxr_linenumber" name="L302" href="#L302">302</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> b = (upper - u) / s;
+<a class="jxr_linenumber" name="L303" href="#L303">303</a> <strong class="jxr_keyword">return</strong> u + moment1(a, b) * s;
+<a class="jxr_linenumber" name="L304" href="#L304">304</a> }
+<a class="jxr_linenumber" name="L305" href="#L305">305</a>
+<a class="jxr_linenumber" name="L306" href="#L306">306</a> <em class="jxr_javadoccomment">/**</em>
+<a class="jxr_linenumber" name="L307" href="#L307">307</a> <em class="jxr_javadoccomment"> * {@inheritDoc}</em>
+<a class="jxr_linenumber" name="L308" href="#L308">308</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L309" href="#L309">309</a> <em class="jxr_javadoccomment"> * <p>Represents the true variance of the truncated normal distribution rather</em>
+<a class="jxr_linenumber" name="L310" href="#L310">310</a> <em class="jxr_javadoccomment"> * than the parent normal distribution variance.</em>
+<a class="jxr_linenumber" name="L311" href="#L311">311</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L312" href="#L312">312</a> <em class="jxr_javadoccomment"> * <p>For \( \mu \) mean of the parent normal distribution,</em>
+<a class="jxr_linenumber" name="L313" href="#L313">313</a> <em class="jxr_javadoccomment"> * \( \sigma \) standard deviation of the parent normal distribution, and</em>
+<a class="jxr_linenumber" name="L314" href="#L314">314</a> <em class="jxr_javadoccomment"> * \( a \lt b \) the truncation interval of the parent normal distribution, the variance is:</em>
+<a class="jxr_linenumber" name="L315" href="#L315">315</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L316" href="#L316">316</a> <em class="jxr_javadoccomment"> * <p>\[ \sigma^2 \left[1 + \frac{a\phi(a)-b\phi(b)}{\Phi(b) - \Phi(a)} -</em>
+<a class="jxr_linenumber" name="L317" href="#L317">317</a> <em class="jxr_javadoccomment"> * \left( \frac{\phi(a)-\phi(b)}{\Phi(b) - \Phi(a)} \right)^2 \right] \]</em>
+<a class="jxr_linenumber" name="L318" href="#L318">318</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L319" href="#L319">319</a> <em class="jxr_javadoccomment"> * <p>where \( \phi \) is the probability density function of the standard normal distribution</em>
+<a class="jxr_linenumber" name="L320" href="#L320">320</a> <em class="jxr_javadoccomment"> * and \( \Phi \) is its cumulative distribution function.</em>
+<a class="jxr_linenumber" name="L321" href="#L321">321</a> <em class="jxr_javadoccomment"> */</em>
+<a class="jxr_linenumber" name="L322" href="#L322">322</a> @Override
+<a class="jxr_linenumber" name="L323" href="#L323">323</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> getVariance() {
+<a class="jxr_linenumber" name="L324" href="#L324">324</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> u = parentNormal.getMean();
+<a class="jxr_linenumber" name="L325" href="#L325">325</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> s = parentNormal.getStandardDeviation();
+<a class="jxr_linenumber" name="L326" href="#L326">326</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> a = (lower - u) / s;
+<a class="jxr_linenumber" name="L327" href="#L327">327</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> b = (upper - u) / s;
+<a class="jxr_linenumber" name="L328" href="#L328">328</a> <strong class="jxr_keyword">return</strong> variance(a, b) * s * s;
+<a class="jxr_linenumber" name="L329" href="#L329">329</a> }
+<a class="jxr_linenumber" name="L330" href="#L330">330</a>
+<a class="jxr_linenumber" name="L331" href="#L331">331</a> <em class="jxr_javadoccomment">/**</em>
+<a class="jxr_linenumber" name="L332" href="#L332">332</a> <em class="jxr_javadoccomment"> * {@inheritDoc}</em>
+<a class="jxr_linenumber" name="L333" href="#L333">333</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L334" href="#L334">334</a> <em class="jxr_javadoccomment"> * <p>The lower bound of the support is equal to the lower bound parameter</em>
+<a class="jxr_linenumber" name="L335" href="#L335">335</a> <em class="jxr_javadoccomment"> * of the distribution.</em>
+<a class="jxr_linenumber" name="L336" href="#L336">336</a> <em class="jxr_javadoccomment"> */</em>
+<a class="jxr_linenumber" name="L337" href="#L337">337</a> @Override
+<a class="jxr_linenumber" name="L338" href="#L338">338</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> getSupportLowerBound() {
+<a class="jxr_linenumber" name="L339" href="#L339">339</a> <strong class="jxr_keyword">return</strong> lower;
+<a class="jxr_linenumber" name="L340" href="#L340">340</a> }
+<a class="jxr_linenumber" name="L341" href="#L341">341</a>
+<a class="jxr_linenumber" name="L342" href="#L342">342</a> <em class="jxr_javadoccomment">/**</em>
+<a class="jxr_linenumber" name="L343" href="#L343">343</a> <em class="jxr_javadoccomment"> * {@inheritDoc}</em>
+<a class="jxr_linenumber" name="L344" href="#L344">344</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L345" href="#L345">345</a> <em class="jxr_javadoccomment"> * <p>The upper bound of the support is equal to the upper bound parameter</em>
+<a class="jxr_linenumber" name="L346" href="#L346">346</a> <em class="jxr_javadoccomment"> * of the distribution.</em>
+<a class="jxr_linenumber" name="L347" href="#L347">347</a> <em class="jxr_javadoccomment"> */</em>
+<a class="jxr_linenumber" name="L348" href="#L348">348</a> @Override
+<a class="jxr_linenumber" name="L349" href="#L349">349</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> getSupportUpperBound() {
+<a class="jxr_linenumber" name="L350" href="#L350">350</a> <strong class="jxr_keyword">return</strong> upper;
+<a class="jxr_linenumber" name="L351" href="#L351">351</a> }
+<a class="jxr_linenumber" name="L352" href="#L352">352</a>
+<a class="jxr_linenumber" name="L353" href="#L353">353</a> <em class="jxr_javadoccomment">/**</em>
+<a class="jxr_linenumber" name="L354" href="#L354">354</a> <em class="jxr_javadoccomment"> * Clip the value to the range [lower, upper].</em>
+<a class="jxr_linenumber" name="L355" href="#L355">355</a> <em class="jxr_javadoccomment"> * This is used to handle floating-point error at the support bound.</em>
+<a class="jxr_linenumber" name="L356" href="#L356">356</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L357" href="#L357">357</a> <em class="jxr_javadoccomment"> * @param x Value x</em>
+<a class="jxr_linenumber" name="L358" href="#L358">358</a> <em class="jxr_javadoccomment"> * @return x clipped to the range</em>
+<a class="jxr_linenumber" name="L359" href="#L359">359</a> <em class="jxr_javadoccomment"> */</em>
+<a class="jxr_linenumber" name="L360" href="#L360">360</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">double</strong> clipToRange(<strong class="jxr_keyword">double</strong> x) {
+<a class="jxr_linenumber" name="L361" href="#L361">361</a> <strong class="jxr_keyword">return</strong> clip(x, lower, upper);
+<a class="jxr_linenumber" name="L362" href="#L362">362</a> }
+<a class="jxr_linenumber" name="L363" href="#L363">363</a>
+<a class="jxr_linenumber" name="L364" href="#L364">364</a> <em class="jxr_javadoccomment">/**</em>
+<a class="jxr_linenumber" name="L365" href="#L365">365</a> <em class="jxr_javadoccomment"> * Clip the value to the range [lower, upper].</em>
+<a class="jxr_linenumber" name="L366" href="#L366">366</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L367" href="#L367">367</a> <em class="jxr_javadoccomment"> * @param x Value x</em>
+<a class="jxr_linenumber" name="L368" href="#L368">368</a> <em class="jxr_javadoccomment"> * @param lower Lower bound (inclusive)</em>
+<a class="jxr_linenumber" name="L369" href="#L369">369</a> <em class="jxr_javadoccomment"> * @param upper Upper bound (inclusive)</em>
+<a class="jxr_linenumber" name="L370" href="#L370">370</a> <em class="jxr_javadoccomment"> * @return x clipped to the range</em>
+<a class="jxr_linenumber" name="L371" href="#L371">371</a> <em class="jxr_javadoccomment"> */</em>
+<a class="jxr_linenumber" name="L372" href="#L372">372</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">double</strong> clip(<strong class="jxr_keyword">double</strong> x, <strong class="jxr_keyword">double</strong> lower, <strong class="jxr_keyword">double</strong> upper) {
+<a class="jxr_linenumber" name="L373" href="#L373">373</a> <strong class="jxr_keyword">if</strong> (x <= lower) {
+<a class="jxr_linenumber" name="L374" href="#L374">374</a> <strong class="jxr_keyword">return</strong> lower;
+<a class="jxr_linenumber" name="L375" href="#L375">375</a> }
+<a class="jxr_linenumber" name="L376" href="#L376">376</a> <strong class="jxr_keyword">return</strong> x < upper ? x : upper;
+<a class="jxr_linenumber" name="L377" href="#L377">377</a> }
+<a class="jxr_linenumber" name="L378" href="#L378">378</a>
+<a class="jxr_linenumber" name="L379" href="#L379">379</a> <em class="jxr_comment">// Calculation of variance and mean can suffer from cancellation.</em>
+<a class="jxr_linenumber" name="L380" href="#L380">380</a> <em class="jxr_comment">//</em>
+<a class="jxr_linenumber" name="L381" href="#L381">381</a> <em class="jxr_comment">// Use formulas from Jorge Fernandez-de-Cossio-Diaz adapted under the</em>
+<a class="jxr_linenumber" name="L382" href="#L382">382</a> <em class="jxr_comment">// terms of the MIT "Expat" License (see NOTICE and LICENSE).</em>
+<a class="jxr_linenumber" name="L383" href="#L383">383</a> <em class="jxr_comment">//</em>
+<a class="jxr_linenumber" name="L384" href="#L384">384</a> <em class="jxr_comment">// These formulas use the complementary error function</em>
+<a class="jxr_linenumber" name="L385" href="#L385">385</a> <em class="jxr_comment">// erfcx(z) = erfc(z) * exp(z^2)</em>
+<a class="jxr_linenumber" name="L386" href="#L386">386</a> <em class="jxr_comment">// This avoids computation of exp terms for the Gaussian PDF and then</em>
+<a class="jxr_linenumber" name="L387" href="#L387">387</a> <em class="jxr_comment">// dividing by the error functions erf or erfc:</em>
+<a class="jxr_linenumber" name="L388" href="#L388">388</a> <em class="jxr_comment">// exp(-0.5*x*x) / erfc(x / sqrt(2)) == 1 / erfcx(x / sqrt(2))</em>
+<a class="jxr_linenumber" name="L389" href="#L389">389</a> <em class="jxr_comment">// At large z the erfcx function is computable but exp(-0.5*z*z) and</em>
+<a class="jxr_linenumber" name="L390" href="#L390">390</a> <em class="jxr_comment">// erfc(z) are zero. Use of these formulas allows computation of the</em>
+<a class="jxr_linenumber" name="L391" href="#L391">391</a> <em class="jxr_comment">// mean and variance for the usable range of the truncated distribution</em>
+<a class="jxr_linenumber" name="L392" href="#L392">392</a> <em class="jxr_comment">// (cdf(a, b) != 0). The variance is not accurate when it approaches</em>
+<a class="jxr_linenumber" name="L393" href="#L393">393</a> <em class="jxr_comment">// machine epsilon (2^-52) at extremely narrow truncations and the</em>
+<a class="jxr_linenumber" name="L394" href="#L394">394</a> <em class="jxr_comment">// computation -> 0.</em>
+<a class="jxr_linenumber" name="L395" href="#L395">395</a> <em class="jxr_comment">//</em>
+<a class="jxr_linenumber" name="L396" href="#L396">396</a> <em class="jxr_comment">// See: https://github.com/cossio/TruncatedNormal.jl</em>
+<a class="jxr_linenumber" name="L397" href="#L397">397</a>
+<a class="jxr_linenumber" name="L398" href="#L398">398</a> <em class="jxr_javadoccomment">/**</em>
+<a class="jxr_linenumber" name="L399" href="#L399">399</a> <em class="jxr_javadoccomment"> * Compute the first moment (mean) of the truncated standard normal distribution.</em>
+<a class="jxr_linenumber" name="L400" href="#L400">400</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L401" href="#L401">401</a> <em class="jxr_javadoccomment"> * <p>Assumes {@code a <= b}.</em>
+<a class="jxr_linenumber" name="L402" href="#L402">402</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L403" href="#L403">403</a> <em class="jxr_javadoccomment"> * @param a Lower bound</em>
+<a class="jxr_linenumber" name="L404" href="#L404">404</a> <em class="jxr_javadoccomment"> * @param b Upper bound</em>
+<a class="jxr_linenumber" name="L405" href="#L405">405</a> <em class="jxr_javadoccomment"> * @return the first moment</em>
+<a class="jxr_linenumber" name="L406" href="#L406">406</a> <em class="jxr_javadoccomment"> */</em>
+<a class="jxr_linenumber" name="L407" href="#L407">407</a> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">double</strong> moment1(<strong class="jxr_keyword">double</strong> a, <strong class="jxr_keyword">double</strong> b) {
+<a class="jxr_linenumber" name="L408" href="#L408">408</a> <em class="jxr_comment">// Assume a <= b</em>
+<a class="jxr_linenumber" name="L409" href="#L409">409</a> <strong class="jxr_keyword">if</strong> (a == b) {
+<a class="jxr_linenumber" name="L410" href="#L410">410</a> <strong class="jxr_keyword">return</strong> a;
+<a class="jxr_linenumber" name="L411" href="#L411">411</a> }
+<a class="jxr_linenumber" name="L412" href="#L412">412</a> <strong class="jxr_keyword">if</strong> (Math.abs(a) > Math.abs(b)) {
+<a class="jxr_linenumber" name="L413" href="#L413">413</a> <em class="jxr_comment">// Subtract from zero to avoid generating -0.0</em>
+<a class="jxr_linenumber" name="L414" href="#L414">414</a> <strong class="jxr_keyword">return</strong> 0 - moment1(-b, -a);
+<a class="jxr_linenumber" name="L415" href="#L415">415</a> }
+<a class="jxr_linenumber" name="L416" href="#L416">416</a>
+<a class="jxr_linenumber" name="L417" href="#L417">417</a> <em class="jxr_comment">// Here:</em>
+<a class="jxr_linenumber" name="L418" href="#L418">418</a> <em class="jxr_comment">// |a| <= |b|</em>
+<a class="jxr_linenumber" name="L419" href="#L419">419</a> <em class="jxr_comment">// a < b</em>
+<a class="jxr_linenumber" name="L420" href="#L420">420</a> <em class="jxr_comment">// 0 < b</em>
+<a class="jxr_linenumber" name="L421" href="#L421">421</a>
+<a class="jxr_linenumber" name="L422" href="#L422">422</a> <strong class="jxr_keyword">if</strong> (a <= -MAX_X) {
+<a class="jxr_linenumber" name="L423" href="#L423">423</a> <em class="jxr_comment">// No truncation</em>
+<a class="jxr_linenumber" name="L424" href="#L424">424</a> <strong class="jxr_keyword">return</strong> 0;
+<a class="jxr_linenumber" name="L425" href="#L425">425</a> }
+<a class="jxr_linenumber" name="L426" href="#L426">426</a> <strong class="jxr_keyword">if</strong> (b >= MAX_X) {
+<a class="jxr_linenumber" name="L427" href="#L427">427</a> <em class="jxr_comment">// One-sided truncation</em>
+<a class="jxr_linenumber" name="L428" href="#L428">428</a> <strong class="jxr_keyword">return</strong> ROOT_2_PI / Erfcx.value(a / ROOT2);
+<a class="jxr_linenumber" name="L429" href="#L429">429</a> }
+<a class="jxr_linenumber" name="L430" href="#L430">430</a>
+<a class="jxr_linenumber" name="L431" href="#L431">431</a> <em class="jxr_comment">// pdf = exp(-0.5*x*x) / sqrt(2*pi)</em>
+<a class="jxr_linenumber" name="L432" href="#L432">432</a> <em class="jxr_comment">// cdf = erfc(-x/sqrt(2)) / 2</em>
+<a class="jxr_linenumber" name="L433" href="#L433">433</a> <em class="jxr_comment">// Compute:</em>
+<a class="jxr_linenumber" name="L434" href="#L434">434</a> <em class="jxr_comment">// -(pdf(b) - pdf(a)) / cdf(b, a)</em>
+<a class="jxr_linenumber" name="L435" href="#L435">435</a> <em class="jxr_comment">// Note:</em>
+<a class="jxr_linenumber" name="L436" href="#L436">436</a> <em class="jxr_comment">// exp(-0.5*b*b) - exp(-0.5*a*a)</em>
+<a class="jxr_linenumber" name="L437" href="#L437">437</a> <em class="jxr_comment">// Use cancellation of powers:</em>
+<a class="jxr_linenumber" name="L438" href="#L438">438</a> <em class="jxr_comment">// exp(-0.5*(b*b-a*a)) * exp(-0.5*a*a) - exp(-0.5*a*a)</em>
+<a class="jxr_linenumber" name="L439" href="#L439">439</a> <em class="jxr_comment">// expm1(-0.5*(b*b-a*a)) * exp(-0.5*a*a)</em>
+<a class="jxr_linenumber" name="L440" href="#L440">440</a>
+<a class="jxr_linenumber" name="L441" href="#L441">441</a> <em class="jxr_comment">// dx = -0.5*(b*b-a*a)</em>
+<a class="jxr_linenumber" name="L442" href="#L442">442</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> dx = 0.5 * (b + a) * (b - a);
+<a class="jxr_linenumber" name="L443" href="#L443">443</a> <strong class="jxr_keyword">double</strong> m;
+<a class="jxr_linenumber" name="L444" href="#L444">444</a> <strong class="jxr_keyword">if</strong> (a <= 0) {
+<a class="jxr_linenumber" name="L445" href="#L445">445</a> <em class="jxr_comment">// Opposite signs</em>
+<a class="jxr_linenumber" name="L446" href="#L446">446</a> m = ROOT_2_PI * -Math.expm1(-dx) * Math.exp(-0.5 * a * a) / ErfDifference.value(a / ROOT2, b / ROOT2);
+<a class="jxr_linenumber" name="L447" href="#L447">447</a> } <strong class="jxr_keyword">else</strong> {
+<a class="jxr_linenumber" name="L448" href="#L448">448</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> z = Math.exp(-dx) * Erfcx.value(b / ROOT2) - Erfcx.value(a / ROOT2);
+<a class="jxr_linenumber" name="L449" href="#L449">449</a> <strong class="jxr_keyword">if</strong> (z == 0) {
+<a class="jxr_linenumber" name="L450" href="#L450">450</a> <em class="jxr_comment">// Occurs when a and b have large magnitudes and are very close</em>
+<a class="jxr_linenumber" name="L451" href="#L451">451</a> <strong class="jxr_keyword">return</strong> (a + b) * 0.5;
+<a class="jxr_linenumber" name="L452" href="#L452">452</a> }
+<a class="jxr_linenumber" name="L453" href="#L453">453</a> m = ROOT_2_PI * Math.expm1(-dx) / z;
+<a class="jxr_linenumber" name="L454" href="#L454">454</a> }
+<a class="jxr_linenumber" name="L455" href="#L455">455</a>
+<a class="jxr_linenumber" name="L456" href="#L456">456</a> <em class="jxr_comment">// Clip to the range</em>
+<a class="jxr_linenumber" name="L457" href="#L457">457</a> <strong class="jxr_keyword">return</strong> clip(m, a, b);
+<a class="jxr_linenumber" name="L458" href="#L458">458</a> }
+<a class="jxr_linenumber" name="L459" href="#L459">459</a>
+<a class="jxr_linenumber" name="L460" href="#L460">460</a> <em class="jxr_javadoccomment">/**</em>
+<a class="jxr_linenumber" name="L461" href="#L461">461</a> <em class="jxr_javadoccomment"> * Compute the second moment of the truncated standard normal distribution.</em>
+<a class="jxr_linenumber" name="L462" href="#L462">462</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L463" href="#L463">463</a> <em class="jxr_javadoccomment"> * <p>Assumes {@code a <= b}.</em>
+<a class="jxr_linenumber" name="L464" href="#L464">464</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L465" href="#L465">465</a> <em class="jxr_javadoccomment"> * @param a Lower bound</em>
+<a class="jxr_linenumber" name="L466" href="#L466">466</a> <em class="jxr_javadoccomment"> * @param b Upper bound</em>
+<a class="jxr_linenumber" name="L467" href="#L467">467</a> <em class="jxr_javadoccomment"> * @return the first moment</em>
+<a class="jxr_linenumber" name="L468" href="#L468">468</a> <em class="jxr_javadoccomment"> */</em>
+<a class="jxr_linenumber" name="L469" href="#L469">469</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">double</strong> moment2(<strong class="jxr_keyword">double</strong> a, <strong class="jxr_keyword">double</strong> b) {
+<a class="jxr_linenumber" name="L470" href="#L470">470</a> <em class="jxr_comment">// Assume a < b.</em>
+<a class="jxr_linenumber" name="L471" href="#L471">471</a> <em class="jxr_comment">// a == b is handled in the variance method</em>
+<a class="jxr_linenumber" name="L472" href="#L472">472</a> <strong class="jxr_keyword">if</strong> (Math.abs(a) > Math.abs(b)) {
+<a class="jxr_linenumber" name="L473" href="#L473">473</a> <strong class="jxr_keyword">return</strong> moment2(-b, -a);
+<a class="jxr_linenumber" name="L474" href="#L474">474</a> }
+<a class="jxr_linenumber" name="L475" href="#L475">475</a>
+<a class="jxr_linenumber" name="L476" href="#L476">476</a> <em class="jxr_comment">// Here:</em>
+<a class="jxr_linenumber" name="L477" href="#L477">477</a> <em class="jxr_comment">// |a| <= |b|</em>
+<a class="jxr_linenumber" name="L478" href="#L478">478</a> <em class="jxr_comment">// a < b</em>
+<a class="jxr_linenumber" name="L479" href="#L479">479</a> <em class="jxr_comment">// 0 < b</em>
+<a class="jxr_linenumber" name="L480" href="#L480">480</a>
+<a class="jxr_linenumber" name="L481" href="#L481">481</a> <strong class="jxr_keyword">if</strong> (a <= -MAX_X) {
+<a class="jxr_linenumber" name="L482" href="#L482">482</a> <em class="jxr_comment">// No truncation</em>
+<a class="jxr_linenumber" name="L483" href="#L483">483</a> <strong class="jxr_keyword">return</strong> 1;
+<a class="jxr_linenumber" name="L484" href="#L484">484</a> }
+<a class="jxr_linenumber" name="L485" href="#L485">485</a> <strong class="jxr_keyword">if</strong> (b >= MAX_X) {
+<a class="jxr_linenumber" name="L486" href="#L486">486</a> <em class="jxr_comment">// One-sided truncation.</em>
+<a class="jxr_linenumber" name="L487" href="#L487">487</a> <em class="jxr_comment">// For a -> inf : moment2 -> a*a</em>
+<a class="jxr_linenumber" name="L488" href="#L488">488</a> <em class="jxr_comment">// This occurs when erfcx(z) is approximated by (1/sqrt(pi)) / z and terms</em>
+<a class="jxr_linenumber" name="L489" href="#L489">489</a> <em class="jxr_comment">// cancel. z > 6.71e7, a > 9.49e7</em>
+<a class="jxr_linenumber" name="L490" href="#L490">490</a> <strong class="jxr_keyword">return</strong> 1 + ROOT_2_PI * a / Erfcx.value(a / ROOT2);
+<a class="jxr_linenumber" name="L491" href="#L491">491</a> }
+<a class="jxr_linenumber" name="L492" href="#L492">492</a>
+<a class="jxr_linenumber" name="L493" href="#L493">493</a> <em class="jxr_comment">// pdf = exp(-0.5*x*x) / sqrt(2*pi)</em>
+<a class="jxr_linenumber" name="L494" href="#L494">494</a> <em class="jxr_comment">// cdf = erfc(-x/sqrt(2)) / 2</em>
+<a class="jxr_linenumber" name="L495" href="#L495">495</a> <em class="jxr_comment">// Compute:</em>
+<a class="jxr_linenumber" name="L496" href="#L496">496</a> <em class="jxr_comment">// 1 - (b*pdf(b) - a*pdf(a)) / cdf(b, a)</em>
+<a class="jxr_linenumber" name="L497" href="#L497">497</a> <em class="jxr_comment">// = (cdf(b, a) - b*pdf(b) -a*pdf(a)) / cdf(b, a)</em>
+<a class="jxr_linenumber" name="L498" href="#L498">498</a>
+<a class="jxr_linenumber" name="L499" href="#L499">499</a> <em class="jxr_comment">// Note:</em>
+<a class="jxr_linenumber" name="L500" href="#L500">500</a> <em class="jxr_comment">// For z -> 0:</em>
+<a class="jxr_linenumber" name="L501" href="#L501">501</a> <em class="jxr_comment">// sqrt(pi / 2) * erf(z / sqrt(2)) -> z</em>
+<a class="jxr_linenumber" name="L502" href="#L502">502</a> <em class="jxr_comment">// z * Math.exp(-0.5 * z * z) -> z</em>
+<a class="jxr_linenumber" name="L503" href="#L503">503</a> <em class="jxr_comment">// Both computations below have cancellation as b -> 0 and the</em>
+<a class="jxr_linenumber" name="L504" href="#L504">504</a> <em class="jxr_comment">// second moment is not computable as the fraction P/Q</em>
+<a class="jxr_linenumber" name="L505" href="#L505">505</a> <em class="jxr_comment">// since P < ulp(Q). This always occurs when b < MIN_X</em>
+<a class="jxr_linenumber" name="L506" href="#L506">506</a> <em class="jxr_comment">// if MIN_X is set at the point where</em>
+<a class="jxr_linenumber" name="L507" href="#L507">507</a> <em class="jxr_comment">// exp(-0.5 * z * z) / sqrt(2 pi) == 1 / sqrt(2 pi).</em>
+<a class="jxr_linenumber" name="L508" href="#L508">508</a> <em class="jxr_comment">// This is JDK dependent due to variations in Math.exp.</em>
+<a class="jxr_linenumber" name="L509" href="#L509">509</a> <em class="jxr_comment">// For b < MIN_X the second moment can be approximated using</em>
+<a class="jxr_linenumber" name="L510" href="#L510">510</a> <em class="jxr_comment">// a uniform distribution: (b^3 - a^3) / (3b - 3a).</em>
+<a class="jxr_linenumber" name="L511" href="#L511">511</a> <em class="jxr_comment">// In practice it also occurs when b > MIN_X since any a < MIN_X</em>
+<a class="jxr_linenumber" name="L512" href="#L512">512</a> <em class="jxr_comment">// is effectively zero for part of the computation. A</em>
+<a class="jxr_linenumber" name="L513" href="#L513">513</a> <em class="jxr_comment">// threshold to transition to a uniform distribution</em>
+<a class="jxr_linenumber" name="L514" href="#L514">514</a> <em class="jxr_comment">// approximation is a compromise. Also note it will not</em>
+<a class="jxr_linenumber" name="L515" href="#L515">515</a> <em class="jxr_comment">// correct computation when (b-a) is small and is far from 0.</em>
+<a class="jxr_linenumber" name="L516" href="#L516">516</a> <em class="jxr_comment">// Thus the second moment is left to be inaccurate for</em>
+<a class="jxr_linenumber" name="L517" href="#L517">517</a> <em class="jxr_comment">// small ranges (b-a) and the variance -> 0 when the true</em>
+<a class="jxr_linenumber" name="L518" href="#L518">518</a> <em class="jxr_comment">// variance is close to or below machine epsilon.</em>
+<a class="jxr_linenumber" name="L519" href="#L519">519</a>
+<a class="jxr_linenumber" name="L520" href="#L520">520</a> <strong class="jxr_keyword">double</strong> m;
+<a class="jxr_linenumber" name="L521" href="#L521">521</a>
+<a class="jxr_linenumber" name="L522" href="#L522">522</a> <strong class="jxr_keyword">if</strong> (a <= 0) {
+<a class="jxr_linenumber" name="L523" href="#L523">523</a> <em class="jxr_comment">// Opposite signs</em>
+<a class="jxr_linenumber" name="L524" href="#L524">524</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> ea = ROOT_PI_2 * Erf.value(a / ROOT2);
+<a class="jxr_linenumber" name="L525" href="#L525">525</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> eb = ROOT_PI_2 * Erf.value(b / ROOT2);
+<a class="jxr_linenumber" name="L526" href="#L526">526</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> fa = ea - a * Math.exp(-0.5 * a * a);
+<a class="jxr_linenumber" name="L527" href="#L527">527</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> fb = eb - b * Math.exp(-0.5 * b * b);
+<a class="jxr_linenumber" name="L528" href="#L528">528</a> <em class="jxr_comment">// Assume fb >= fa && eb >= ea</em>
+<a class="jxr_linenumber" name="L529" href="#L529">529</a> <em class="jxr_comment">// If fb <= fa this is a tiny range around 0</em>
+<a class="jxr_linenumber" name="L530" href="#L530">530</a> m = (fb - fa) / (eb - ea);
+<a class="jxr_linenumber" name="L531" href="#L531">531</a> <em class="jxr_comment">// Clip to the range</em>
+<a class="jxr_linenumber" name="L532" href="#L532">532</a> m = clip(m, 0, 1);
+<a class="jxr_linenumber" name="L533" href="#L533">533</a> } <strong class="jxr_keyword">else</strong> {
+<a class="jxr_linenumber" name="L534" href="#L534">534</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> dx = 0.5 * (b + a) * (b - a);
+<a class="jxr_linenumber" name="L535" href="#L535">535</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> ex = Math.exp(-dx);
+<a class="jxr_linenumber" name="L536" href="#L536">536</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> ea = ROOT_PI_2 * Erfcx.value(a / ROOT2);
+<a class="jxr_linenumber" name="L537" href="#L537">537</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> eb = ROOT_PI_2 * Erfcx.value(b / ROOT2);
+<a class="jxr_linenumber" name="L538" href="#L538">538</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> fa = ea + a;
+<a class="jxr_linenumber" name="L539" href="#L539">539</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> fb = eb + b;
+<a class="jxr_linenumber" name="L540" href="#L540">540</a> m = (fa - fb * ex) / (ea - eb * ex);
+<a class="jxr_linenumber" name="L541" href="#L541">541</a> <em class="jxr_comment">// Clip to the range</em>
+<a class="jxr_linenumber" name="L542" href="#L542">542</a> m = clip(m, a * a, b * b);
+<a class="jxr_linenumber" name="L543" href="#L543">543</a> }
+<a class="jxr_linenumber" name="L544" href="#L544">544</a> <strong class="jxr_keyword">return</strong> m;
+<a class="jxr_linenumber" name="L545" href="#L545">545</a> }
+<a class="jxr_linenumber" name="L546" href="#L546">546</a>
+<a class="jxr_linenumber" name="L547" href="#L547">547</a> <em class="jxr_javadoccomment">/**</em>
+<a class="jxr_linenumber" name="L548" href="#L548">548</a> <em class="jxr_javadoccomment"> * Compute the variance of the truncated standard normal distribution.</em>
+<a class="jxr_linenumber" name="L549" href="#L549">549</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L550" href="#L550">550</a> <em class="jxr_javadoccomment"> * <p>Assumes {@code a <= b}.</em>
+<a class="jxr_linenumber" name="L551" href="#L551">551</a> <em class="jxr_javadoccomment"> *</em>
+<a class="jxr_linenumber" name="L552" href="#L552">552</a> <em class="jxr_javadoccomment"> * @param a Lower bound</em>
+<a class="jxr_linenumber" name="L553" href="#L553">553</a> <em class="jxr_javadoccomment"> * @param b Upper bound</em>
+<a class="jxr_linenumber" name="L554" href="#L554">554</a> <em class="jxr_javadoccomment"> * @return the first moment</em>
+<a class="jxr_linenumber" name="L555" href="#L555">555</a> <em class="jxr_javadoccomment"> */</em>
+<a class="jxr_linenumber" name="L556" href="#L556">556</a> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">double</strong> variance(<strong class="jxr_keyword">double</strong> a, <strong class="jxr_keyword">double</strong> b) {
+<a class="jxr_linenumber" name="L557" href="#L557">557</a> <strong class="jxr_keyword">if</strong> (a == b) {
+<a class="jxr_linenumber" name="L558" href="#L558">558</a> <strong class="jxr_keyword">return</strong> 0;
+<a class="jxr_linenumber" name="L559" href="#L559">559</a> }
+<a class="jxr_linenumber" name="L560" href="#L560">560</a>
+<a class="jxr_linenumber" name="L561" href="#L561">561</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> m1 = moment1(a, b);
+<a class="jxr_linenumber" name="L562" href="#L562">562</a> <strong class="jxr_keyword">double</strong> m2 = moment2(a, b);
+<a class="jxr_linenumber" name="L563" href="#L563">563</a> <em class="jxr_comment">// variance = m2 - m1*m1</em>
+<a class="jxr_linenumber" name="L564" href="#L564">564</a> <em class="jxr_comment">// rearrange x^2 - y^2 as (x-y)(x+y)</em>
+<a class="jxr_linenumber" name="L565" href="#L565">565</a> m2 = Math.sqrt(m2);
+<a class="jxr_linenumber" name="L566" href="#L566">566</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> variance = (m2 - m1) * (m2 + m1);
+<a class="jxr_linenumber" name="L567" href="#L567">567</a>
+<a class="jxr_linenumber" name="L568" href="#L568">568</a> <em class="jxr_comment">// Detect floating-point error.</em>
+<a class="jxr_linenumber" name="L569" href="#L569">569</a> <strong class="jxr_keyword">if</strong> (variance >= 1) {
+<a class="jxr_linenumber" name="L570" href="#L570">570</a> <em class="jxr_comment">// Note:</em>
+<a class="jxr_linenumber" name="L571" href="#L571">571</a> <em class="jxr_comment">// Extreme truncations in the tails can compute a variance above 1,</em>
+<a class="jxr_linenumber" name="L572" href="#L572">572</a> <em class="jxr_comment">// for example if m2 is infinite: m2 - m1*m1 > 1</em>
+<a class="jxr_linenumber" name="L573" href="#L573">573</a> <em class="jxr_comment">// Detect no truncation as the terms a and b lie far either side of zero;</em>
+<a class="jxr_linenumber" name="L574" href="#L574">574</a> <em class="jxr_comment">// otherwise return 0 to indicate very small unknown variance.</em>
+<a class="jxr_linenumber" name="L575" href="#L575">575</a> <strong class="jxr_keyword">return</strong> a < -1 && b > 1 ? 1 : 0;
+<a class="jxr_linenumber" name="L576" href="#L576">576</a> } <strong class="jxr_keyword">else</strong> <strong class="jxr_keyword">if</strong> (variance <= 0) {
+<a class="jxr_linenumber" name="L577" href="#L577">577</a> <em class="jxr_comment">// Floating-point error can create negative variance so return 0.</em>
+<a class="jxr_linenumber" name="L578" href="#L578">578</a> <strong class="jxr_keyword">return</strong> 0;
+<a class="jxr_linenumber" name="L579" href="#L579">579</a> }
+<a class="jxr_linenumber" name="L580" href="#L580">580</a>
+<a class="jxr_linenumber" name="L581" href="#L581">581</a> <strong class="jxr_keyword">return</strong> variance;
+<a class="jxr_linenumber" name="L582" href="#L582">582</a> }
+<a class="jxr_linenumber" name="L583" href="#L583">583</a> }
+</pre>
+<hr/>
+<div id="footer">Copyright © 2018–2022 <a href="https://www.apache.org/">The Apache Software Foundation</a>. All rights reserved.</div>
+</body>
+</html>