001 /*
002 * Java Genetic Algorithm Library (jenetics-4.3.0).
003 * Copyright (c) 2007-2018 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 3.1
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 getCount() {
114 return _n;
115 }
116
117 /**
118 * Return the arithmetic mean of values recorded, or {@code Double.NaN} if
119 * no values have been recorded.
120 *
121 * @return the arithmetic mean of values, or zero if none
122 */
123 public double getMean() {
124 return _n == 0L ? NaN : _m1.value();
125 }
126
127 /**
128 * Return the variance of values recorded, or {@code Double.NaN} if no
129 * values have been recorded.
130 *
131 * @return the variance of values, or {@code NaN} if none
132 */
133 public double getVariance() {
134 double var = NaN;
135 if (_n == 1L) {
136 var = _m2.value();
137 } else if (_n > 1L) {
138 var = _m2.value()/(_n - 1.0);
139 }
140
141 return var;
142 }
143
144 /**
145 * Return the skewness of values recorded, or {@code Double.NaN} if less
146 * than two values have been recorded.
147 *
148 * @see <a href="https://en.wikipedia.org/wiki/Skewness">Skewness</a>
149 *
150 * @return the skewness of values, or {@code NaN} if less than two values
151 * have been recorded
152 */
153 public double getSkewness() {
154 double skewness = NaN;
155 if (_n >= 3L) {
156 final double var = _m2.value()/(_n - 1.0);
157 if (var < 10E-20) {
158 skewness = 0.0d;
159 } else {
160 skewness = (_n*_m3.value())/
161 ((_n - 1.0)*(_n - 2.0)*sqrt(var)*var);
162 }
163 }
164
165 return skewness;
166 }
167
168 /**
169 * Return the kurtosis of values recorded, or {@code Double.NaN} if less
170 * than four values have been recorded.
171 *
172 * @see <a href="https://en.wikipedia.org/wiki/Kurtosis">Kurtosis</a>
173 *
174 * @return the kurtosis of values, or {@code NaN} if less than four values
175 * have been recorded
176 */
177 public double getKurtosis() {
178 double kurtosis = NaN;
179 if (_n > 3L) {
180 final double var = _m2.value()/(_n - 1);
181 if (_n <= 3L || var < 10E-20) {
182 kurtosis = 0.0;
183 } else {
184 kurtosis = (_n*(_n + 1.0)*_m4.value() -
185 3.0*_m2.value()*_m2.value()*(_n - 1.0))/
186 ((_n - 1.0)*(_n - 2.0)*(_n - 3.0)*var*var);
187 }
188 }
189 return kurtosis;
190 }
191
192 final boolean sameState(final MomentStatistics statistics) {
193 return _n == statistics._n &&
194 _m1.sameState(statistics._m1) &&
195 _m2.sameState(statistics._m2) &&
196 _m3.sameState(statistics._m3) &&
197 _m4.sameState(statistics._m4);
198 }
199
200 }
|