1 /******************************************************************************* 2 3 Struct implementation of a moving average designed to handle irregularly 4 spaced data. In principle this can also be used as a serializable record 5 struct. 6 7 Copyright: 8 Copyright (c) 2009-2016 dunnhumby Germany GmbH. 9 All rights reserved. 10 11 License: 12 Boost Software License Version 1.0. See LICENSE_BOOST.txt for details. 13 Alternatively, this file may be distributed under the terms of the Tango 14 3-Clause BSD License (see LICENSE_BSD.txt for details). 15 16 *******************************************************************************/ 17 18 module ocean.math.IrregularMovingAverage; 19 20 21 22 import ocean.core.Enforce; 23 24 import ocean.meta.traits.Basic /* isIntegerType, isFloatingPointType */; 25 import ocean.math.IEEE; 26 import ocean.math.Math; 27 import core.stdc.time; 28 29 import ocean.core.Verify; 30 31 version (unittest) import ocean.core.Test; 32 33 /******************************************************************************* 34 35 Moving average designed to handle irregularly spaced data, i.e. data where 36 the time intervals between successive values are not the same. 37 38 Params: 39 Time = type used to indicate the observation time of individual 40 data values; may be any type implicitly convertible to 41 time_t, or any integral or floating point value 42 Value = floating-point type used to indicate data values 43 44 The double-exponential moving average implemented here is derived from 45 section 3 (eqs. 16-22) of Cipra, T. (2006) 'Exponential smoothing for 46 irregular data', Applications of Mathematics 51 (6): 597-604 47 <http://dml.cz/handle/10338.dmlcz/134655>. 48 49 Most of the parameter names are derived from their counterparts in the 50 aforementioned article. 51 52 *******************************************************************************/ 53 54 public struct IrregularMovingAverage (Time = time_t, Value = real) 55 { 56 static assert(is(Time : time_t) 57 || isIntegerType!(Time) 58 || isFloatingPointType!(Time)); 59 60 static assert(isFloatingPointType!(Value)); 61 62 63 /*************************************************************************** 64 65 Internal helper alias 66 67 ***************************************************************************/ 68 69 private alias typeof(this) This; 70 71 72 /*************************************************************************** 73 74 Time this moving average was last updated 75 76 ***************************************************************************/ 77 78 private Time update_time_; 79 80 81 /*************************************************************************** 82 83 Discount coefficient that determines how much weight will be given to 84 new data values. Must be in the range 0 < beta < 1 (note the strict 85 inequality). 86 87 ***************************************************************************/ 88 89 private Value beta_; 90 91 92 /*************************************************************************** 93 94 Time-dependent smoothing coefficient 95 96 ***************************************************************************/ 97 98 private Value alpha; 99 100 101 /*************************************************************************** 102 103 Auxiliary values used to help calculate the moving average 104 105 ***************************************************************************/ 106 107 private Value w; 108 /// ditto 109 private Value z; 110 111 112 /*************************************************************************** 113 114 Single- and double-exponential smoothing statistics 115 116 ***************************************************************************/ 117 118 private Value s; 119 /// ditto 120 private Value s2; 121 122 123 /*************************************************************************** 124 125 (In)sanity checks on internal data 126 127 ***************************************************************************/ 128 129 invariant () 130 { 131 assert(!isNaN(this.beta_)); 132 assert(!isInfinity(this.beta_)); 133 assert(0 < this.beta_); 134 assert(this.beta_ < 1); 135 136 assert(!isNaN(this.alpha)); 137 assert(!isInfinity(this.alpha)); 138 assert(0 <= this.alpha); 139 assert(this.alpha <= 1); 140 141 assert(!isNaN(this.w)); 142 assert(!isInfinity(this.w)); 143 144 assert(!isNaN(this.z)); 145 assert(!isInfinity(this.z)); 146 147 assert(!isNaN(this.s)); 148 assert(!isInfinity(this.s)); 149 150 assert(!isNaN(this.s2)); 151 assert(!isInfinity(this.s2)); 152 153 static if (isFloatingPointType!(Time)) 154 { 155 assert(!isNaN(this.update_time_)); 156 assert(!isInfinity(this.update_time_)); 157 } 158 } 159 160 161 162 /*************************************************************************** 163 164 "Constructor"-style static opCall 165 166 Params: 167 beta = discount coefficient that determines how much contribution 168 new values make to the average (must have 0 < beta < 1) 169 expected_time_interval = expected average time interval between 170 successive data points (must be > 0) 171 first_value = initial value of the quantity being averaged 172 first_value_time = observation time of the value provided 173 174 Returns: 175 IrregularMovingAverage instance based on the single provided 176 data point 177 178 ***************************************************************************/ 179 180 public static This opCall (Value beta, Time expected_time_interval, 181 Value first_value, Time first_value_time) 182 out (ima) 183 { 184 assert(&ima); // confirm the invariant holds 185 } 186 do 187 { 188 verify(!isNaN(beta)); 189 verify(!isInfinity(beta)); 190 verify(0.0 < beta); 191 verify(beta < 1.0); 192 193 verify(!isNaN(first_value)); 194 verify(!isInfinity(first_value)); 195 196 static if (isFloatingPointType!(Time)) 197 { 198 verify(!isNaN(expected_time_interval)); 199 verify(!isInfinity(expected_time_interval)); 200 201 verify(!isNaN(first_value_time)); 202 verify(!isInfinity(first_value_time)); 203 } 204 205 verify(expected_time_interval > 0); 206 207 // cast to real is necessary to satisfy 'pow' overrides 208 Value beta_pow = pow(beta, cast(real) expected_time_interval); 209 Value alpha_start = 1.0 - beta_pow; 210 Value wz_start = 211 (alpha_start * alpha_start) / (expected_time_interval * beta_pow); 212 213 This ima = 214 { 215 update_time_: first_value_time, 216 beta_: beta, 217 alpha: alpha_start, 218 w: wz_start, 219 z: wz_start, 220 s: first_value, 221 s2: first_value 222 }; 223 224 return ima; 225 } 226 227 228 /*************************************************************************** 229 230 Generates an updated moving average with a new data point, using the 231 specified discount coefficient. The new moving average is returned 232 as a new struct instance, rather than mutating the internal struct 233 data. 234 235 Params: 236 value = new data value to include in the average 237 value_time = observation time of the value provided 238 239 ***************************************************************************/ 240 241 public void update (Value value, Time value_time) 242 { 243 verify(!isNaN(value)); 244 verify(!isInfinity(value)); 245 246 static if (isFloatingPointType!(Time)) 247 { 248 verify(!isNaN(value_time)); 249 verify(!isInfinity(value_time)); 250 } 251 252 enforce(value_time > this.update_time, 253 "Cannot incorporate data point from the past!"); 254 255 // values re-used multiple times updating parameters 256 Value time_diff = value_time - this.update_time; 257 Value beta_pow = pow(this.beta, time_diff); 258 259 this.update_time_ = value_time; 260 261 // update w based using old alpha value 262 this.w = this.w / 263 (beta_pow + (time_diff * beta_pow * this.w / this.alpha)); 264 265 this.alpha = this.alpha / (beta_pow + this.alpha); 266 267 // update z using new alpha and w values 268 this.z = this.z / 269 (beta_pow + (this.alpha * this.z / this.w)); 270 271 // update s using new alpha value 272 this.s = (this.alpha * value) + ((1.0 - this.alpha) * this.s); 273 274 // update s2 using new alpha and new s value 275 this.s2 = (this.alpha * this.s) + ((1.0 - this.alpha) * this.s2); 276 } 277 278 279 /*************************************************************************** 280 281 Returns: 282 the current double-exponentially-smoothed moving average 283 284 ***************************************************************************/ 285 286 public Value value () 287 { 288 return this.s2; 289 } 290 291 292 /*************************************************************************** 293 294 Calculates the predicted future value of the moving average 295 296 Params: 297 t = time for which to calculate the predicted value; must be 298 greater than or equal to the update time of the average 299 (if not provided, current time will be used) 300 301 Returns: 302 the predicted value of the averaged quantity at time t 303 304 Throws: 305 exception if t is less than the update time of the average 306 307 ***************************************************************************/ 308 309 public Value predicted_value (Time t) 310 out (estimate) 311 { 312 assert(!isNaN(estimate)); 313 assert(!isInfinity(estimate)); 314 } 315 do 316 { 317 enforce(t >= this.update_time_, 318 "Cannot calculate predicted value from the past!"); 319 320 Time k = t - this.update_time; 321 322 assert(k >= 0); 323 324 Value estimate = this.s + 325 ( ((this.z / this.w) + ((this.z / this.alpha) * k)) * (this.s - this.s2) ); 326 327 return estimate; 328 } 329 330 331 /*************************************************************************** 332 333 Returns: 334 most recent update time of the moving average (i.e. the time 335 of the last data point included) 336 337 ***************************************************************************/ 338 339 public Time update_time () 340 { 341 return this.update_time_; 342 } 343 344 345 /*************************************************************************** 346 347 Private setter for update_time, used internally only 348 349 Params: 350 t = new update time to set; must be greater than existing value 351 352 Returns: 353 the newly-set update time value (i.e. t) 354 355 ***************************************************************************/ 356 357 private Time update_time (Time t) 358 { 359 static if (isFloatingPointType!(Time)) 360 { 361 verify(!isNaN(value_time)); 362 verify(!isInfinity(value_time)); 363 } 364 365 verify(t > this.update_time_); 366 367 return this.update_time_ = t; 368 } 369 370 371 /*************************************************************************** 372 373 Returns: 374 currently-set discount coefficient beta for this moving average 375 376 ***************************************************************************/ 377 378 public Value beta () 379 { 380 return this.beta_; 381 } 382 383 384 /*************************************************************************** 385 386 Set the discount coefficient for this moving average. It may be 387 useful to adjust this value if, for example, the average waiting 388 time between data points turns out to be very different from 389 initial expectations. 390 391 Params: 392 b = new value to set for the discount coefficient beta; must 393 have 0 < b < 1. 394 395 Returns: 396 newly-set discount coefficient beta for this moving average 397 (i.e. b). 398 399 ***************************************************************************/ 400 401 public Value beta (Value b) 402 { 403 verify(!isNaN(b)); 404 verify(!isInfinity(b)); 405 verify(0 < b); 406 verify(b < 1); 407 408 return this.beta_ = b; 409 } 410 } 411 412 unittest 413 { 414 // --- basic initialization tests --- 415 416 real val0 = 23.5; 417 time_t t0 = 123456; 418 419 auto avg0 = IrregularMovingAverage!()(0.5, 60, val0, t0); 420 421 test(avg0.value is val0); 422 test(avg0.update_time == t0); 423 424 /* Without multiple statistics, cannot make 425 * predictions about the future 426 */ 427 test(avg0.predicted_value(t0) == val0); 428 test(avg0.predicted_value(time(null)) == val0); 429 test(avg0.predicted_value(123456789) == val0); 430 431 432 // --- test updates where discount coefficient is altered --- 433 434 real val1 = 36.5; 435 time_t t1 = 123459; 436 437 /* Create a copy and update with an 438 * increased discount coefficient 439 */ 440 auto avg1a = avg0; 441 avg1a.beta = 0.9; 442 avg1a.update(val1, t1); 443 444 // larger data second value ==> increase in the average 445 test(avg1a.update_time == t1); 446 test(avg1a.value > val0); 447 test(avg1a.value < val1); 448 449 time_t t_future = 1234560; 450 451 /* Predicted value for update time should 452 * be equal to the observed value, but note 453 * that precise equality here is extremely 454 * unlikely due to floating point arithmetic. 455 * For this reason, these tests are disabled 456 * until we identify a nice "approximately 457 * equal" method to use. 458 */ 459 //assert(avg1a.predicted_value(t1) is val1); 460 //assert(avg1a.predicted_value(avg1a.update_time) is val1); 461 462 // should now have an increasing trend 463 test(avg1a.predicted_value(t_future) > avg1a.value); 464 test(avg1a.predicted_value(t_future) > avg1a.predicted_value(t1)); 465 466 /* Now we create another copy of avg0 and update 467 * with a _decreased_ discount coefficient, but 468 * the same second data point as avg1a 469 */ 470 auto avg1b = avg0; 471 avg1b.beta = 0.1; 472 avg1b.update(val1, t1); 473 474 // larger second data value ==> increase in the average 475 test(avg1b.update_time == t1); 476 test(avg1b.value > val0); 477 test(avg1b.value < val1); 478 479 /* Predicted value for update time should 480 * be equal to the observed value, but note 481 * that precise equality here is extremely 482 * unlikely due to floating point arithmetic. 483 * For this reason, these tests are disabled 484 * until we identify a nice "approximately 485 * equal" method to use. 486 */ 487 //assert(avg1b.predicted_value(t1) is val1); 488 //assert(avg1b.predicted_value(avg1b.update_time) is val1); 489 490 // should now have an increasing trend 491 test(avg1b.predicted_value(t_future) > avg1b.value); 492 test(avg1b.predicted_value(t_future) > avg1b.predicted_value(t1)); 493 494 /* greater weight to newer records should 495 * make for larger values 496 */ 497 test(avg1b.value > avg1a.value); 498 499 500 // --- test updates where discount coefficient is not altered --- 501 502 real val2 = 19.5; 503 time_t t2 = t1; 504 505 /* Create a copy of avg0 and update with 506 * a new (smaller-valued) data point 507 */ 508 auto avg2a = avg0; 509 avg2a.update(val2, t2); 510 511 // smaller second value ==> decrease in the average 512 test(avg2a.update_time == t2); 513 test(avg2a.value < val0); 514 test(avg2a.value > val2); 515 516 // should now have an decreasing trend 517 test(avg2a.predicted_value(t_future) < avg2a.value); 518 test(avg2a.predicted_value(t_future) < avg2a.predicted_value(t2)); 519 520 /* should get identical results if we repeat 521 * the same moving-average update 522 */ 523 auto avg2b = avg0; 524 avg2b.update(val2, t2); 525 test(avg2a == avg2b); 526 }