001 /*
002 * Java Genetic Algorithm Library (jenetics-5.2.0).
003 * Copyright (c) 2007-2020 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@gmail.com)
019 */
020 package io.jenetics.stat;
021
022 import static java.lang.Double.NaN;
023 import static java.lang.Math.sqrt;
024 import static java.util.Objects.requireNonNull;
025
026 import io.jenetics.internal.math.DoubleAdder;
027
028 /**
029 * Base class for statistical moments calculation.
030 *
031 * @see <a href="http://people.xiph.org/~tterribe/notes/homs.html">
032 * Computing Higher-Order Moments Online</a>
033 *
034 * @author <a href="mailto:franz.wilhelmstoetter@gmail.com">Franz Wilhelmstötter</a>
035 * @since 3.0
036 * @version 5.2
037 */
038 abstract class MomentStatistics {
039
040 // the sample count.
041 private long _n = 0L;
042
043 // Variables used for statistical moments.
044 private final DoubleAdder _m1 = new DoubleAdder();
045 private final DoubleAdder _m2 = new DoubleAdder();
046 private final DoubleAdder _m3 = new DoubleAdder();
047 private final DoubleAdder _m4 = new DoubleAdder();
048
049
050 /**
051 * Update the moments with the given {@code value}.
052 *
053 * @param value the value which is used to update this statistical moments.
054 */
055 void accept(final double value) {
056 ++_n;
057
058 final double n = _n;
059 final double d = value - _m1.value();
060 final double dN = d/n;
061 final double dN2 = dN*dN;
062 final double t1 = d*dN*(n - 1.0);
063
064 _m1.add(dN);
065 _m4.add(t1*dN2*(n*n - 3.0*n + 3.0))
066 .add(6.0*dN2*_m2.value() - 4.0*dN*_m3.value());
067 _m3.add(t1*dN*(n - 2.0) - 3.0*dN*_m2.value());
068 _m2.add(t1);
069 }
070
071 /**
072 * Combines the state of another {@code Moments} object into this one.
073 *
074 * @see <a href="http://people.xiph.org/~tterribe/notes/homs.html">
075 * Computing Higher-Order Moments Online</a>
076 */
077 void combine(final MomentStatistics b) {
078 requireNonNull(b);
079
080 final double m2 = _m2.value();
081 final double m3 = _m3.value();
082
083 final double pn = _n;
084 final double n = _n + b._n;
085 final double nn = n*n;
086
087 final double d = b._m1.value() - _m1.value();
088 final double dd = d*d;
089
090 _n += b._n;
091
092 _m1.add(d*b._n/n);
093
094 _m2.add(b._m2).add(dd*pn*b._n/n);
095
096 _m3.add(b._m3)
097 .add(dd*d*(pn*b._n*(pn - b._n)/nn))
098 .add(3.0*d*(pn*b._m2.value() - b._n*m2)/n);
099
100 final double bnbn = (double)b._n*(double)b._n;
101
102 _m4.add(b._m4)
103 .add(dd*dd*(pn*b._n*(pn*pn - pn*b._n + bnbn)/(nn*n)))
104 .add(6.0*dd*(pn*pn*b._m2.value() + bnbn*m2)/nn)
105 .add(4.0*d*(pn*b._m3.value() - b._n*m3)/n);
106 }
107
108 /**
109 * Returns the count of values recorded.
110 *
111 * @return the count of recorded values
112 */
113 public long count() {
114 return _n;
115 }
116
117 /**
118 * Returns the count of values recorded.
119 *
120 * @return the count of recorded values
121 * @deprecated Use {@link #count()} instead
122 */
123 @Deprecated
124 public long getCount() {
125 return _n;
126 }
127
128 /**
129 * Return the arithmetic mean of values recorded, or {@code Double.NaN} if
130 * no values have been recorded.
131 *
132 * @return the arithmetic mean of values, or zero if none
133 */
134 public double mean() {
135 return _n == 0L ? NaN : _m1.value();
136 }
137
138 /**
139 * Return the arithmetic mean of values recorded, or {@code Double.NaN} if
140 * no values have been recorded.
141 *
142 * @return the arithmetic mean of values, or zero if none
143 * @deprecated Use {@link #mean()} instead
144 */
145 @Deprecated
146 public double getMean() {
147 return mean();
148 }
149
150 /**
151 * Return the variance of values recorded, or {@code Double.NaN} if no
152 * values have been recorded.
153 *
154 * @return the variance of values, or {@code NaN} if none
155 */
156 public double variance() {
157 double var = NaN;
158 if (_n == 1L) {
159 var = _m2.value();
160 } else if (_n > 1L) {
161 var = _m2.value()/(_n - 1.0);
162 }
163
164 return var;
165 }
166
167 /**
168 * Return the variance of values recorded, or {@code Double.NaN} if no
169 * values have been recorded.
170 *
171 * @return the variance of values, or {@code NaN} if none
172 * @deprecated Use {@link #variance()} instead
173 */
174 @Deprecated
175 public double getVariance() {
176 return variance();
177 }
178
179 /**
180 * Return the skewness of values recorded, or {@code Double.NaN} if less
181 * than two values have been recorded.
182 *
183 * @see <a href="https://en.wikipedia.org/wiki/Skewness">Skewness</a>
184 *
185 * @return the skewness of values, or {@code NaN} if less than two values
186 * have been recorded
187 */
188 public double skewness() {
189 double skewness = NaN;
190 if (_n >= 3L) {
191 final double var = _m2.value()/(_n - 1.0);
192 skewness = var < 10E-20
193 ? 0.0d
194 : (_n*_m3.value())/((_n - 1.0)*(_n - 2.0)*sqrt(var)*var);
195 }
196
197 return skewness;
198 }
199
200 /**
201 * Return the skewness of values recorded, or {@code Double.NaN} if less
202 * than two values have been recorded.
203 *
204 * @see <a href="https://en.wikipedia.org/wiki/Skewness">Skewness</a>
205 *
206 * @return the skewness of values, or {@code NaN} if less than two values
207 * have been recorded
208 * @deprecated Use {@link #skewness()} instead
209 */
210 @Deprecated
211 public double getSkewness() {
212 return skewness();
213 }
214
215 /**
216 * Return the kurtosis of values recorded, or {@code Double.NaN} if less
217 * than four values have been recorded.
218 *
219 * @see <a href="https://en.wikipedia.org/wiki/Kurtosis">Kurtosis</a>
220 *
221 * @return the kurtosis of values, or {@code NaN} if less than four values
222 * have been recorded
223 */
224 public double kurtosis() {
225 double kurtosis = NaN;
226 if (_n > 3L) {
227 final double var = _m2.value()/(_n - 1);
228 kurtosis = _n <= 3L || var < 10E-20
229 ? 0.0
230 : (_n*(_n + 1.0)*_m4.value() -
231 3.0*_m2.value()*_m2.value()*(_n - 1.0))/
232 ((_n - 1.0)*(_n - 2.0)*(_n - 3.0)*var*var);
233 }
234 return kurtosis;
235 }
236
237 /**
238 * Return the kurtosis of values recorded, or {@code Double.NaN} if less
239 * than four values have been recorded.
240 *
241 * @see <a href="https://en.wikipedia.org/wiki/Kurtosis">Kurtosis</a>
242 *
243 * @return the kurtosis of values, or {@code NaN} if less than four values
244 * have been recorded
245 * @deprecated Use {@link #kurtosis()} instead
246 */
247 @Deprecated
248 public double getKurtosis() {
249 return kurtosis();
250 }
251
252 final boolean sameState(final MomentStatistics statistics) {
253 return _n == statistics._n &&
254 _m1.sameState(statistics._m1) &&
255 _m2.sameState(statistics._m2) &&
256 _m3.sameState(statistics._m3) &&
257 _m4.sameState(statistics._m4);
258 }
259
260 }
|