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