001 /*
002 * Java Genetic Algorithm Library (jenetics-2.0.2).
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 2.0 — <em>$Date: 2014-03-28 $</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 public double getValue() {
117 return _q[2];
118 }
119
120 @Override
121 public void accumulate(final N value) {
122 if (!_initialized) {
123 initialize(value.doubleValue());
124 } else {
125 update(value.doubleValue());
126 }
127
128 ++_samples;
129 }
130
131
132 private void initialize(double value) {
133 if (_n[0] < 0.0) {
134 _n[0] = 0.0;
135 _q[0] = value;
136 } else if (_n[1] == 0.0) {
137 _n[1] = 1.0;
138 _q[1] = value;
139 } else if (_n[2] == 0.0) {
140 _n[2] = 2.0;
141 _q[2] = value;
142 } else if (_n[3] == 0.0) {
143 _n[3] = 3.0;
144 _q[3] = value;
145 } else if (_n[4] == 0.0) {
146 _n[4] = 4.0;
147 _q[4] = value;
148 }
149
150 if (_n[4] != 0.0) {
151 Arrays.sort(_q);
152
153 _nn[0] = 2.0*_quantile;
154 _nn[1] = 4.0*_quantile;
155 _nn[2] = 2.0*_quantile + 2.0;
156
157 _dn[0] = _quantile/2.0;
158 _dn[1] = _quantile;
159 _dn[2] = (1.0 + _quantile)/2.0;
160
161 _initialized = true;
162 }
163 }
164
165 private void update(double value) {
166 assert (_initialized);
167
168 // If min or max, handle as special case; otherwise, ...
169 if (_quantile == 0.0) {
170 if (value < _q[2]) {
171 _q[2] = value;
172 }
173 } else if (_quantile == 1.0) {
174 if (value > _q[2]) {
175 _q[2] = value;
176 }
177 } else {
178 // Increment marker locations and update min and max.
179 if (value < _q[0]) {
180 ++_n[1]; ++_n[2]; ++_n[3]; ++_n[4]; _q[0] = value;
181 } else if (value < _q[1]) {
182 ++_n[1]; ++_n[2]; ++_n[3]; ++_n[4];
183 } else if (value < _q[2]) {
184 ++_n[2]; ++_n[3]; ++_n[4];
185 } else if (value < _q[3]) {
186 ++_n[3]; ++_n[4];
187 } else if (value < _q[4]) {
188 ++_n[4];
189 } else {
190 ++_n[4]; _q[4] = value;
191 }
192
193 // Increment positions of markers k + 1
194 _nn[0] += _dn[0];
195 _nn[1] += _dn[1];
196 _nn[2] += _dn[2];
197
198 // Adjust heights of markers 0 to 2 if necessary
199 double mm = _n[1] - 1.0;
200 double mp = _n[1] + 1.0;
201 if (_nn[0] >= mp && _n[2] > mp) {
202 _q[1] = qPlus(mp, _n[0], _n[1], _n[2], _q[0], _q[1], _q[2]);
203 _n[1] = mp;
204 } else if (_nn[0] <= mm && _n[0] < mm) {
205 _q[1] = qMinus(mm, _n[0], _n[1], _n[2], _q[0], _q[1], _q[2]);
206 _n[1] = mm;
207 }
208
209 mm = _n[2] - 1.0;
210 mp = _n[2] + 1.0;
211 if (_nn[1] >= mp && _n[3] > mp) {
212 _q[2] = qPlus(mp, _n[1], _n[2], _n[3], _q[1], _q[2], _q[3]);
213 _n[2] = mp;
214 } else if (_nn[1] <= mm && _n[1] < mm) {
215 _q[2] = qMinus(mm, _n[1], _n[2], _n[3], _q[1], _q[2], _q[3]);
216 _n[2] = mm;
217 }
218
219 mm = _n[3] - 1.0;
220 mp = _n[3] + 1.0;
221 if (_nn[2] >= mp && _n[4] > mp) {
222 _q[3] = qPlus(mp, _n[2], _n[3], _n[4], _q[2], _q[3], _q[4]);
223 _n[3] = mp;
224 } else if (_nn[2] <= mm && _n[2] < mm) {
225 _q[3] = qMinus(mm, _n[2], _n[3], _n[4], _q[2], _q[3], _q[4]);
226 _n[3] = mm;
227 }
228 }
229 }
230
231 private static double qPlus(
232 final double mp,
233 final double m0,
234 final double m1,
235 final double m2,
236 final double q0,
237 final double q1,
238 final double q2
239 ) {
240 double result = q1 +
241 ((mp - m0)*(q2 - q1)/(m2 - m1) +
242 (m2 - mp)*(q1 - q0)/(m1 - m0))/(m2 - m0);
243
244 if (result > q2) {
245 result = q1 + (q2 - q1)/(m2 - m1);
246 }
247
248 return result;
249 }
250
251 private static double qMinus(
252 final double mm,
253 final double m0,
254 final double m1,
255 final double m2,
256 final double q0,
257 final double q1,
258 final double q2
259 ) {
260 double result = q1 -
261 ((mm - m0)*(q2 - q1)/(m2 - m1) +
262 (m2 - mm)*(q1 - q0)/(m1 - m0))/(m2 - m0);
263
264 if (q0 > result) {
265 result = q1 + (q0 - q1)/(m0 - m1);
266 }
267
268 return result;
269 }
270
271 @Override
272 public int hashCode() {
273 return HashBuilder.of(getClass()).
274 and(super.hashCode()).
275 and(_quantile).
276 and(_dn).
277 and(_n).
278 and(_nn).
279 and(_q).value();
280 }
281
282 @Override
283 public boolean equals(final Object obj) {
284 if (obj == this) {
285 return true;
286 }
287 if (obj == null || getClass() != obj.getClass()) {
288 return false;
289 }
290
291 final Quantile<?> quantile = (Quantile<?>)obj;
292 return 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 @Override
301 public String toString() {
302 return format(
303 "%s[samples=%d, quantile=%f]",
304 getClass().getSimpleName(), getSamples(), getValue()
305 );
306 }
307
308 @Override
309 public Quantile<N> clone() {
310 return (Quantile<N>)super.clone();
311 }
312
313
314 static <N extends Number> Quantile<N> median() {
315 return new Quantile<>(0.5);
316 }
317
318 }
|