1 /*
2 * Copyright (C) 2015 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.statistics;
17
18 /**
19 * Defines the error function and methods related to it.
20 * The error function (or Gaussian function) is a special function typically
21 * used in statistics and probability or to solve partial differential
22 * equations.
23 * The error function is defined as the integral between a given interval of an
24 * exponential with a negative squared exponent (like gaussian functions).
25 * Such integral can be used for purposes like obtaining the probability of
26 * gaussian distributions, because such probability is the integral of a pdf
27 * having a similar expression as the error function.
28 * The error function (a.k.a. erf) has the following properties:
29 * - is an odd function erf(-x) = -erf(x)
30 * - erf(0) = 0
31 * - erf(infinity) = 1
32 * <p>
33 * Typically, also the complementary error function (a.k.a. erfc) is also used,
34 * and it is defined as erfc(x) = 1 - erf(x)
35 * The complementary error function has the following properties:
36 * - erfc(0) = 1
37 * - erfc(infinity) = 0
38 * - erfc(-x) = 2 - erfc(x).
39 * <p>
40 * Both functions (erf and erfc) are special cases of the incomplete gamma
41 * function.
42 * <p>
43 * This class is based in code of Numerical Recipes 3rd ed. section 6.2.2.
44 */
45 public class Erf {
46 /**
47 * Number of coefficients to approximate the error function using a
48 * Chebychev approximation.
49 */
50 private static final int N_COF = 28;
51
52 /**
53 * Coefficients of Chebychev polynomial to approximate the error function.
54 */
55 private static final double[] COF = {-1.3026537197817094,
56 6.4196979235649026e-1, 1.9476473204185836e-2, -9.561514786808631e-3,
57 -9.46595344482036e-4, 3.66839497852761e-4, 4.2523324806907e-5,
58 -2.0278578112534e-5, -1.624290004647e-6, 1.303655835580e-6,
59 1.5626441722e-8, -8.5238095915e-8, 6.529054439e-9, 5.059343495e-9,
60 -9.91364156e-10, -2.27365122e-10, 9.6467911e-11, 2.394038e-12,
61 -6.886027e-12, 8.94487e-13, 3.13092e-13, -1.12708e-13, 3.81e-16,
62 7.106e-15, -1.523e-15, -9.4e-17, 1.21e-16, -2.8e-17};
63
64 /**
65 * Empty constructor.
66 */
67 private Erf() {
68 }
69
70 /**
71 * Evaluates the error function at x.
72 *
73 * @param x value to evaluate the error function at.
74 * @return value of the error function.
75 */
76 public static double erf(final double x) {
77 if (x >= 0.0) {
78 return 1.0 - erfccheb(x);
79 } else {
80 return erfccheb(-x) - 1.0;
81 }
82 }
83
84 /**
85 * Evaluates the complementary error function at x.
86 *
87 * @param x value to evaluate the complementary error function at.
88 * @return value of the complementary error function.
89 */
90 public static double erfc(final double x) {
91 if (x >= 0.0) {
92 return erfccheb(x);
93 } else {
94 return 2.0 - erfccheb(-x);
95 }
96 }
97
98 /**
99 * Evaluates the inverse of the complementary error function at p.
100 * Then:
101 * p = erfc(x) <--> x = inverfc(p)
102 *
103 * @param p value to evaluate the inverse erfc function at.
104 * @return result of the evaluation.
105 */
106 public static double inverfc(final double p) {
107 double x;
108 double err;
109 final double t;
110 final double pp;
111 if (p >= 2.0) {
112 return -100.0;
113 }
114 if (p <= 0.0) {
115 return 100.0;
116 }
117 pp = p < 1.0 ? p : 2.0 - p;
118 t = Math.sqrt(-2.0 * Math.log(pp / 2.0));
119 x = -0.70711 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t);
120 for (int j = 0; j < 2; j++) {
121 err = erfc(x) - pp;
122 x += err / (1.12837916709551257 * Math.exp(-x * x) - x * err);
123 }
124 return (p < 1.0 ? x : -x);
125 }
126
127 /**
128 * Evaluates the inverse of the error function at p.
129 * Then:
130 * p = erf(x) <--> x = inverf(p)
131 *
132 * @param p value to evaluate the inverse erf function at.
133 * @return result of the evaluation.
134 */
135 public static double inverf(final double p) {
136 return inverfc(1.0 - p);
137 }
138
139 /**
140 * Computes the complementary error function by using the Chebychev method
141 * approximation.
142 *
143 * @param z value to evaluate the function at.
144 * @return evaluation of the erfc at provided value.
145 * @throws IllegalArgumentException if provided value is negative.
146 */
147 private static double erfccheb(final double z) {
148 int j;
149 final double t;
150 final double ty;
151 double tmp;
152 double d = 0.0;
153 double dd = 0.0;
154 if (z < 0.0) {
155 throw new IllegalArgumentException("erfccheb requires non-negative argument");
156 }
157
158 t = 2.0 / (2.0 + z);
159 ty = 4.0 * t - 2.;
160 for (j = N_COF - 1; j > 0; j--) {
161 tmp = d;
162 d = ty * d - dd + COF[j];
163 dd = tmp;
164 }
165 return t * Math.exp(-z * z + 0.5 * (COF[0] + ty * d) - dd);
166 }
167 }