001 /*
002 * Java Genetic Algorithm Library (jenetics-3.0.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.stat;
021
022 import static java.lang.Double.compare;
023 import static java.lang.String.format;
024 import static org.jenetics.internal.util.Equality.eq;
025
026 import java.util.Arrays;
027 import java.util.function.DoubleConsumer;
028
029 import org.jenetics.internal.util.Equality;
030 import org.jenetics.internal.util.Hash;
031
032
033 /**
034 * Implementation of the quantile estimation algorithm published by
035 * <p>
036 * <strong>Raj JAIN and Imrich CHLAMTAC</strong>:
037 * <em>
038 * The P<sup>2</sup> Algorithm for Dynamic Calculation of Quantiles and
039 * Histograms Without Storing Observations
040 * </em>
041 * <br>
042 * [<a href="http://www.cse.wustl.edu/~jain/papers/ftp/psqr.pdf">Communications
043 * of the ACM; October 1985, Volume 28, Number 10</a>]
044 * <p>
045 * <strong>Note that this implementation is not synchronized.</strong> If
046 * multiple threads access this object concurrently, and at least one of the
047 * threads modifies it, it must be synchronized externally.
048 *
049 * @see <a href="http://en.wikipedia.org/wiki/Quantile">Wikipedia: Quantile</a>
050 *
051 * @author <a href="mailto:franz.wilhelmstoetter@gmx.at">Franz Wilhelmstötter</a>
052 * @since 1.0
053 * @version 3.0 — <em>$Date: 2014-12-28 $</em>
054 */
055 public class Quantile implements DoubleConsumer {
056
057 private long _samples = 0;
058
059 // The desired quantile.
060 private final double _quantile;
061
062 // Marker heights.
063 private final double[] _q = {0, 0, 0, 0, 0};
064
065 // Marker positions.
066 private final double[] _n = {0, 0, 0, 0, 0};
067
068 // Desired marker positions.
069 private final double[] _nn = {0, 0, 0};
070
071 // Desired marker position increments.
072 private final double[] _dn = {0, 0, 0};
073
074 private boolean _initialized;
075
076 /**
077 * Create a new quantile accumulator with the given value.
078 *
079 * @param quantile the wished quantile value.
080 * @throws IllegalArgumentException if the {@code quantile} is not in the
081 * range {@code [0, 1]}.
082 */
083 public Quantile(final double quantile) {
084 _quantile = quantile;
085 init(quantile);
086 }
087
088 private void init(final double quantile) {
089 if (quantile < 0.0 || quantile > 1) {
090 throw new IllegalArgumentException(format(
091 "Quantile (%s) not in the valid range of [0, 1]", quantile
092 ));
093 }
094
095 Arrays.fill(_q, 0);
096 Arrays.fill(_n, 0);
097 Arrays.fill(_nn, 0);
098 Arrays.fill(_dn, 0);
099
100 _n[0] = -1.0;
101 _q[2] = 0.0;
102 _initialized = compare(quantile, 0.0) == 0 ||
103 compare(quantile, 1.0) == 0;
104 _samples = 0;
105 }
106
107 /**
108 * Reset this object to its initial state.
109 */
110 public void reset() {
111 init(_quantile);
112 }
113
114 /**
115 * Return the computed quantile value.
116 *
117 * @return the quantile value.
118 */
119 public double getValue() {
120 return _q[2];
121 }
122
123 public long getSamples() {
124 return _samples;
125 }
126
127 @Override
128 public void accept(final double value) {
129 if (!_initialized) {
130 initialize(value);
131 } else {
132 update(value);
133 }
134
135 ++_samples;
136 }
137
138
139 private void initialize(double value) {
140 if (_n[0] < 0.0) {
141 _n[0] = 0.0;
142 _q[0] = value;
143 } else if (_n[1] == 0.0) {
144 _n[1] = 1.0;
145 _q[1] = value;
146 } else if (_n[2] == 0.0) {
147 _n[2] = 2.0;
148 _q[2] = value;
149 } else if (_n[3] == 0.0) {
150 _n[3] = 3.0;
151 _q[3] = value;
152 } else if (_n[4] == 0.0) {
153 _n[4] = 4.0;
154 _q[4] = value;
155 }
156
157 if (_n[4] != 0.0) {
158 Arrays.sort(_q);
159
160 _nn[0] = 2.0*_quantile;
161 _nn[1] = 4.0*_quantile;
162 _nn[2] = 2.0*_quantile + 2.0;
163
164 _dn[0] = _quantile/2.0;
165 _dn[1] = _quantile;
166 _dn[2] = (1.0 + _quantile)/2.0;
167
168 _initialized = true;
169 }
170 }
171
172 private void update(double value) {
173 assert (_initialized);
174
175 // If min or max, handle as special case; otherwise, ...
176 if (_quantile == 0.0) {
177 if (value < _q[2]) {
178 _q[2] = value;
179 }
180 } else if (_quantile == 1.0) {
181 if (value > _q[2]) {
182 _q[2] = value;
183 }
184 } else {
185 // Increment marker locations and update min and max.
186 if (value < _q[0]) {
187 ++_n[1]; ++_n[2]; ++_n[3]; ++_n[4]; _q[0] = value;
188 } else if (value < _q[1]) {
189 ++_n[1]; ++_n[2]; ++_n[3]; ++_n[4];
190 } else if (value < _q[2]) {
191 ++_n[2]; ++_n[3]; ++_n[4];
192 } else if (value < _q[3]) {
193 ++_n[3]; ++_n[4];
194 } else if (value < _q[4]) {
195 ++_n[4];
196 } else {
197 ++_n[4]; _q[4] = value;
198 }
199
200 // Increment positions of markers k + 1
201 _nn[0] += _dn[0];
202 _nn[1] += _dn[1];
203 _nn[2] += _dn[2];
204
205 // Adjust heights of markers 0 to 2 if necessary
206 double mm = _n[1] - 1.0;
207 double mp = _n[1] + 1.0;
208 if (_nn[0] >= mp && _n[2] > mp) {
209 _q[1] = qPlus(mp, _n[0], _n[1], _n[2], _q[0], _q[1], _q[2]);
210 _n[1] = mp;
211 } else if (_nn[0] <= mm && _n[0] < mm) {
212 _q[1] = qMinus(mm, _n[0], _n[1], _n[2], _q[0], _q[1], _q[2]);
213 _n[1] = mm;
214 }
215
216 mm = _n[2] - 1.0;
217 mp = _n[2] + 1.0;
218 if (_nn[1] >= mp && _n[3] > mp) {
219 _q[2] = qPlus(mp, _n[1], _n[2], _n[3], _q[1], _q[2], _q[3]);
220 _n[2] = mp;
221 } else if (_nn[1] <= mm && _n[1] < mm) {
222 _q[2] = qMinus(mm, _n[1], _n[2], _n[3], _q[1], _q[2], _q[3]);
223 _n[2] = mm;
224 }
225
226 mm = _n[3] - 1.0;
227 mp = _n[3] + 1.0;
228 if (_nn[2] >= mp && _n[4] > mp) {
229 _q[3] = qPlus(mp, _n[2], _n[3], _n[4], _q[2], _q[3], _q[4]);
230 _n[3] = mp;
231 } else if (_nn[2] <= mm && _n[2] < mm) {
232 _q[3] = qMinus(mm, _n[2], _n[3], _n[4], _q[2], _q[3], _q[4]);
233 _n[3] = mm;
234 }
235 }
236 }
237
238 private static double qPlus(
239 final double mp,
240 final double m0,
241 final double m1,
242 final double m2,
243 final double q0,
244 final double q1,
245 final double q2
246 ) {
247 double result = q1 +
248 ((mp - m0)*(q2 - q1)/(m2 - m1) +
249 (m2 - mp)*(q1 - q0)/(m1 - m0))/(m2 - m0);
250
251 if (result > q2) {
252 result = q1 + (q2 - q1)/(m2 - m1);
253 }
254
255 return result;
256 }
257
258 private static double qMinus(
259 final double mm,
260 final double m0,
261 final double m1,
262 final double m2,
263 final double q0,
264 final double q1,
265 final double q2
266 ) {
267 double result = q1 -
268 ((mm - m0)*(q2 - q1)/(m2 - m1) +
269 (m2 - mm)*(q1 - q0)/(m1 - m0))/(m2 - m0);
270
271 if (q0 > result) {
272 result = q1 + (q0 - q1)/(m0 - m1);
273 }
274
275 return result;
276 }
277
278 @Override
279 public int hashCode() {
280 return Hash.of(getClass()).
281 and(super.hashCode()).
282 and(_quantile).
283 and(_dn).
284 and(_n).
285 and(_nn).
286 and(_q).value();
287 }
288
289 @Override
290 public boolean equals(final Object obj) {
291 return Equality.of(this, obj).test(quantile ->
292 super.equals(obj) &&
293 eq(_quantile, quantile._quantile) &&
294 eq(_dn, quantile._dn) &&
295 eq(_n, quantile._n) &&
296 eq(_nn, quantile._nn) &&
297 eq(_q, quantile._q)
298 );
299 }
300
301 @Override
302 public String toString() {
303 return format(
304 "%s[samples=%d, quantile=%f]",
305 getClass().getSimpleName(), getSamples(), getValue()
306 );
307 }
308
309
310 static Quantile median() {
311 return new Quantile(0.5);
312 }
313
314 }
|