001 /*
002 * Java Genetic Algorithm Library (jenetics-1.6.0).
003 * Copyright (c) 2007-2014 Franz Wilhelmstötter
004 *
005 * Licensed under the Apache License, Version 2.0 (the "License");
006 * you may not use this file except in compliance with the License.
007 * You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 *
017 * Author:
018 * Franz Wilhelmstötter (franz.wilhelmstoetter@gmx.at)
019 */
020 package org.jenetics.internal.math;
021
022 import static java.lang.Math.PI;
023 import static java.lang.Math.abs;
024 import static java.lang.Math.exp;
025 import static java.lang.Math.log;
026 import static java.lang.Math.sqrt;
027
028 import org.jenetics.util.StaticObject;
029
030 /**
031 * Some statistical special functions.
032 *
033 * @author <a href="mailto:franz.wilhelmstoetter@gmx.at">Franz Wilhelmstötter</a>
034 * @since 1.0
035 * @version 1.4 — <em>$Date: 2014-02-15 $</em>
036 */
037 public final class statistics extends StaticObject {
038 private statistics() {}
039
040
041 /**
042 * Uses Lanczos approximation formula. See Numerical Recipes 6.1.
043 *
044 * @param x the value to calculate the log gamma for.
045 * @return the log gamma value.
046 */
047 public static double logGamma(final double x) {
048 final double ser = 1.0 + 76.18009173 /
049 (x + 0) - 86.50532033 /
050 (x + 1) + 24.01409822 /
051 (x + 2) - 1.231739516 /
052 (x + 3) + 0.00120858003 /
053 (x + 4) - 0.00000536382 /
054 (x + 5);
055
056 return (x - 0.5)*log(x + 4.5) - (x + 4.5) + log(ser*sqrt(2 * PI));
057 }
058
059 public static double gamma(final double x) {
060 return exp(logGamma(x));
061 }
062
063 public static double Γ(final double x) {
064 return gamma(x);
065 }
066
067
068
069 /**
070 * Return the <i>error function</i> of {@code z}. The fractional error of
071 * this implementation is less than 1.2E-7.
072 *
073 * @param z the value to calculate the error function for.
074 * @return the error function for {@code z}.
075 */
076 public static double erf(final double z) {
077 final double t = 1.0/(1.0 + 0.5*abs(z));
078
079 // Horner's method
080 final double result = 1 - t*exp(
081 -z*z - 1.26551223 +
082 t*( 1.00002368 +
083 t*( 0.37409196 +
084 t*( 0.09678418 +
085 t*(-0.18628806 +
086 t*( 0.27886807 +
087 t*(-1.13520398 +
088 t*( 1.48851587 +
089 t*(-0.82215223 +
090 t*(0.17087277))))))))));
091
092 return z >= 0 ? result : -result;
093 }
094
095 /**
096 * Return φ(x), the standard Gaussian pdf.
097 *
098 * @see #φ(double)
099 * @param x the value to calculate φ for.
100 * @return the φ value for x.
101 */
102 public static double phi(final double x) {
103 return exp(-x*x/2.0) / sqrt(2.0*PI);
104 }
105
106 /**
107 * Return φ(x), the standard Gaussian pdf.
108 *
109 * @see #phi(double)
110 * @param x the value to calculate φ for.
111 * @return the φ value for x.
112 */
113 public static double φ(final double x) {
114 return phi(x);
115 }
116
117 /**
118 * Return φ(x, µ, σ), the standard Gaussian pdf with mean µ and stddev σ.
119 *
120 * @see #phi(double, double, double)
121 * @param x the value to calculate φ for.
122 * @param mu the mean value.
123 * @param sigma the stddev.
124 * @return the φ value for x.
125 */
126 public static double phi(final double x, final double mu, final double sigma) {
127 return phi((x - mu)/sigma)/sigma;
128 }
129
130 /**
131 * Return φ(x, µ, σ), the standard Gaussian pdf with mean µ and stddev σ.
132 *
133 * @see #phi(double, double, double)
134 * @param x the value to calculate φ for.
135 * @param µ the mean value.
136 * @param σ the stddev.
137 * @return the φ value for x.
138 */
139 public static double φ(final double x, final double µ, final double σ) {
140 return phi(x, µ, σ);
141 }
142
143 /**
144 * Return Φ(z), the standard Gaussian cdf using Taylor approximation.
145 *
146 * @param z the value to calculate Φ for.
147 * @return the Φ for value z.
148 */
149 public static double Phi(final double z) {
150 if (z < -8.0) {
151 return 0.0;
152 }
153 if (z > 8.0) {
154 return 1.0;
155 }
156
157 double s = 0.0;
158 double t = z;
159 for (int i = 3; s + t != s; i += 2) {
160 s = s + t;
161 t = t*z*z/i;
162 }
163 return 0.5 + s*phi(z);
164 }
165
166 /**
167 * Return Φ(z), the standard Gaussian cdf using Taylor approximation.
168 *
169 * @param z the value to calculate Φ for.
170 * @return the Φ for value z.
171 */
172 public static double Φ(final double z) {
173 return Phi(z);
174 }
175
176 /**
177 * Return Φ(z, µ, σ), the standard Gaussian cdf with mean µ and stddev σ.
178 *
179 * @see #phi(double, double, double)
180 * @param z the value to calculate Φ for.
181 * @param mu the mean value.
182 * @param sigma the stddev.
183 * @return the φ value for x.
184 */
185 public static double Phi(final double z, final double mu, final double sigma) {
186 return Phi((z - mu)/sigma);
187 }
188
189 /**
190 * Return Φ(z, µ, σ), the standard Gaussian cdf with mean µ and stddev σ.
191 *
192 * @see #phi(double, double, double)
193 * @param z the value to calculate Φ for.
194 * @param µ the mean value.
195 * @param σ the stddev.
196 * @return the φ value for x.
197 */
198 public static double Φ(final double z, final double µ, final double σ) {
199 return Phi(z, µ, σ);
200 }
201
202 }
|