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