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.stat;
021
022 import static java.lang.Double.compare;
023 import static java.lang.String.format;
024 import static org.jenetics.internal.util.object.eq;
025
026 import java.util.Arrays;
027
028 import org.jenetics.internal.util.HashBuilder;
029
030 import org.jenetics.util.MappedAccumulator;
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 1.3 — <em>$Date: 2014-03-01 $</em>
054 */
055 public class Quantile<N extends Number> extends MappedAccumulator<N> {
056
057 // The desired quantile.
058 private final double _quantile;
059
060 // Marker heights.
061 private final double[] _q = {0, 0, 0, 0, 0};
062
063 // Marker positions.
064 private final double[] _n = {0, 0, 0, 0, 0};
065
066 // Desired marker positions.
067 private final double[] _nn = {0, 0, 0};
068
069 // Desired marker position increments.
070 private final double[] _dn = {0, 0, 0};
071
072 private boolean _initialized;
073
074 /**
075 * Create a new quantile accumulator with the given value.
076 *
077 * @param quantile the wished quantile value.
078 * @throws IllegalArgumentException if the {@code quantile} is not in the
079 * range {@code [0, 1]}.
080 */
081 public Quantile(final double quantile) {
082 _quantile = quantile;
083 init(quantile);
084 }
085
086 private void init(final double quantile) {
087 if (quantile < 0.0 || quantile > 1) {
088 throw new IllegalArgumentException(format(
089 "Quantile (%s) not in the valid range of [0, 1]", quantile
090 ));
091 }
092
093 Arrays.fill(_q, 0);
094 Arrays.fill(_n, 0);
095 Arrays.fill(_nn, 0);
096 Arrays.fill(_dn, 0);
097
098 _n[0] = -1.0;
099 _q[2] = 0.0;
100 _initialized = compare(quantile, 0.0) == 0 || compare(quantile, 1.0) == 0;
101 _samples = 0;
102 }
103
104 /**
105 * Reset this object to its initial state.
106 */
107 public void reset() {
108 init(_quantile);
109 }
110
111 /**
112 * Return the computed quantile value.
113 *
114 * @return the quantile value.
115 *
116 * @deprecated Use {@link #getValue()} instead.
117 */
118 @Deprecated
119 public double getQuantile() {
120 return _q[2];
121 }
122
123 /**
124 * Return the computed quantile value.
125 *
126 * @return the quantile value.
127 */
128 public double getValue() {
129 return _q[2];
130 }
131
132 @Override
133 public void accumulate(final N value) {
134 if (!_initialized) {
135 initialize(value.doubleValue());
136 } else {
137 update(value.doubleValue());
138 }
139
140 ++_samples;
141 }
142
143
144 private void initialize(double value) {
145 if (_n[0] < 0.0) {
146 _n[0] = 0.0;
147 _q[0] = value;
148 } else if (_n[1] == 0.0) {
149 _n[1] = 1.0;
150 _q[1] = value;
151 } else if (_n[2] == 0.0) {
152 _n[2] = 2.0;
153 _q[2] = value;
154 } else if (_n[3] == 0.0) {
155 _n[3] = 3.0;
156 _q[3] = value;
157 } else if (_n[4] == 0.0) {
158 _n[4] = 4.0;
159 _q[4] = value;
160 }
161
162 if (_n[4] != 0.0) {
163 Arrays.sort(_q);
164
165 _nn[0] = 2.0*_quantile;
166 _nn[1] = 4.0*_quantile;
167 _nn[2] = 2.0*_quantile + 2.0;
168
169 _dn[0] = _quantile/2.0;
170 _dn[1] = _quantile;
171 _dn[2] = (1.0 + _quantile)/2.0;
172
173 _initialized = true;
174 }
175 }
176
177 private void update(double value) {
178 assert (_initialized);
179
180 // If min or max, handle as special case; otherwise, ...
181 if (_quantile == 0.0) {
182 if (value < _q[2]) {
183 _q[2] = value;
184 }
185 } else if (_quantile == 1.0) {
186 if (value > _q[2]) {
187 _q[2] = value;
188 }
189 } else {
190 // Increment marker locations and update min and max.
191 if (value < _q[0]) {
192 ++_n[1]; ++_n[2]; ++_n[3]; ++_n[4]; _q[0] = value;
193 } else if (value < _q[1]) {
194 ++_n[1]; ++_n[2]; ++_n[3]; ++_n[4];
195 } else if (value < _q[2]) {
196 ++_n[2]; ++_n[3]; ++_n[4];
197 } else if (value < _q[3]) {
198 ++_n[3]; ++_n[4];
199 } else if (value < _q[4]) {
200 ++_n[4];
201 } else {
202 ++_n[4]; _q[4] = value;
203 }
204
205 // Increment positions of markers k + 1
206 _nn[0] += _dn[0];
207 _nn[1] += _dn[1];
208 _nn[2] += _dn[2];
209
210 // Adjust heights of markers 0 to 2 if necessary
211 double mm = _n[1] - 1.0;
212 double mp = _n[1] + 1.0;
213 if (_nn[0] >= mp && _n[2] > mp) {
214 _q[1] = qPlus(mp, _n[0], _n[1], _n[2], _q[0], _q[1], _q[2]);
215 _n[1] = mp;
216 } else if (_nn[0] <= mm && _n[0] < mm) {
217 _q[1] = qMinus(mm, _n[0], _n[1], _n[2], _q[0], _q[1], _q[2]);
218 _n[1] = mm;
219 }
220
221 mm = _n[2] - 1.0;
222 mp = _n[2] + 1.0;
223 if (_nn[1] >= mp && _n[3] > mp) {
224 _q[2] = qPlus(mp, _n[1], _n[2], _n[3], _q[1], _q[2], _q[3]);
225 _n[2] = mp;
226 } else if (_nn[1] <= mm && _n[1] < mm) {
227 _q[2] = qMinus(mm, _n[1], _n[2], _n[3], _q[1], _q[2], _q[3]);
228 _n[2] = mm;
229 }
230
231 mm = _n[3] - 1.0;
232 mp = _n[3] + 1.0;
233 if (_nn[2] >= mp && _n[4] > mp) {
234 _q[3] = qPlus(mp, _n[2], _n[3], _n[4], _q[2], _q[3], _q[4]);
235 _n[3] = mp;
236 } else if (_nn[2] <= mm && _n[2] < mm) {
237 _q[3] = qMinus(mm, _n[2], _n[3], _n[4], _q[2], _q[3], _q[4]);
238 _n[3] = mm;
239 }
240 }
241 }
242
243 private static double qPlus(
244 final double mp,
245 final double m0,
246 final double m1,
247 final double m2,
248 final double q0,
249 final double q1,
250 final double q2
251 ) {
252 double result = q1 +
253 ((mp - m0)*(q2 - q1)/(m2 - m1) +
254 (m2 - mp)*(q1 - q0)/(m1 - m0))/(m2 - m0);
255
256 if (result > q2) {
257 result = q1 + (q2 - q1)/(m2 - m1);
258 }
259
260 return result;
261 }
262
263 private static double qMinus(
264 final double mm,
265 final double m0,
266 final double m1,
267 final double m2,
268 final double q0,
269 final double q1,
270 final double q2
271 ) {
272 double result = q1 -
273 ((mm - m0)*(q2 - q1)/(m2 - m1) +
274 (m2 - mm)*(q1 - q0)/(m1 - m0))/(m2 - m0);
275
276 if (q0 > result) {
277 result = q1 + (q0 - q1)/(m0 - m1);
278 }
279
280 return result;
281 }
282
283 @Override
284 public int hashCode() {
285 return HashBuilder.of(getClass()).
286 and(super.hashCode()).
287 and(_quantile).
288 and(_dn).
289 and(_n).
290 and(_nn).
291 and(_q).value();
292 }
293
294 @Override
295 public boolean equals(final Object obj) {
296 if (obj == this) {
297 return true;
298 }
299 if (obj == null || getClass() != obj.getClass()) {
300 return false;
301 }
302
303 final Quantile<?> quantile = (Quantile<?>)obj;
304 return super.equals(obj) &&
305 eq(_quantile, quantile._quantile) &&
306 eq(_dn, quantile._dn) &&
307 eq(_n, quantile._n) &&
308 eq(_nn, quantile._nn) &&
309 eq(_q, quantile._q);
310 }
311
312 @Override
313 public String toString() {
314 return format(
315 "%s[samples=%d, quantile=%f]",
316 getClass().getSimpleName(), getSamples(), getQuantile()
317 );
318 }
319
320 @Override
321 public Quantile<N> clone() {
322 return (Quantile<N>)super.clone();
323 }
324
325
326 static <N extends Number> Quantile<N> median() {
327 return new Quantile<>(0.5);
328 }
329
330 }
|