timeseries.cpp
Go to the documentation of this file.
00001 /*************************************************************************** 00002 file : $URL: https://frepple.svn.sourceforge.net/svnroot/frepple/trunk/modules/forecast/forecastsolver.cpp $ 00003 version : $LastChangedRevision: 1505 $ $LastChangedBy: jdetaeye $ 00004 date : $LastChangedDate: 2011-08-26 18:55:08 +0200 (Fri, 26 Aug 2011) $ 00005 ***************************************************************************/ 00006 00007 /*************************************************************************** 00008 * * 00009 * Copyright (C) 2007-2011 by Johan De Taeye, frePPLe bvba * 00010 * * 00011 * This library is free software; you can redistribute it and/or modify it * 00012 * under the terms of the GNU Lesser General Public License as published * 00013 * by the Free Software Foundation; either version 2.1 of the License, or * 00014 * (at your option) any later version. * 00015 * * 00016 * This library is distributed in the hope that it will be useful, * 00017 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 00018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser * 00019 * General Public License for more details. * 00020 * * 00021 * You should have received a copy of the GNU Lesser General Public * 00022 * License along with this library; if not, write to the Free Software * 00023 * Foundation Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,* 00024 * USA * 00025 * * 00026 ***************************************************************************/ 00027 00028 #include "forecast.h" 00029 00030 namespace module_forecast 00031 { 00032 00033 #define ACCURACY 0.01 00034 00035 void Forecast::generateFutureValues( 00036 const double history[], unsigned int historycount, 00037 const Date buckets[], unsigned int bucketcount, 00038 bool debug) 00039 { 00040 // Validate the input 00041 if (!history || !buckets) 00042 throw RuntimeException("Null argument to forecast function"); 00043 if (bucketcount < 2) 00044 throw DataException("Need at least 2 forecast dates"); 00045 00046 // Strip zero demand buckets at the start. 00047 // Eg when demand only starts in the middle of horizon, we only want to 00048 // use the second part of the horizon for forecasting. The zeros before the 00049 // demand start would distort the results. 00050 while (historycount >= 1 && *history == 0.0) 00051 { 00052 ++history; 00053 --historycount; 00054 } 00055 00056 // We create the forecasting objects in stack memory for best performance. 00057 MovingAverage moving_avg; 00058 Croston croston; 00059 SingleExponential single_exp; 00060 DoubleExponential double_exp; 00061 Seasonal seasonal; 00062 int numberOfMethods = 4; 00063 ForecastMethod* methods[4]; 00064 00065 // Rules to determine which forecast methods can be applied 00066 methods[0] = &moving_avg; 00067 if (historycount < getForecastSkip() + 5) 00068 // Too little history: only use moving average 00069 numberOfMethods = 1; 00070 else 00071 { 00072 unsigned int zero = 0; 00073 for (unsigned long i = 0; i < historycount; ++i) 00074 if (!history[i]) ++zero; 00075 if (zero > Croston::getMinIntermittence() * historycount) 00076 { 00077 // If there are too many zeros: use croston or moving average. 00078 numberOfMethods = 2; 00079 methods[1] = &croston; 00080 } 00081 else 00082 { 00083 // The most common case: enough values and not intermittent 00084 methods[1] = &single_exp; 00085 methods[2] = &double_exp; 00086 methods[3] = &seasonal; 00087 } 00088 } 00089 00090 // Initialize a vector with the smape weights 00091 double *weight = new double[historycount+1]; 00092 weight[historycount] = 1.0; 00093 for (int i=historycount-1; i>=0; --i) 00094 weight[i] = weight[i+1] * Forecast::getForecastSmapeAlfa(); 00095 00096 // Evaluate each forecast method 00097 double best_error = DBL_MAX; 00098 int best_method = -1; 00099 double error; 00100 try 00101 { 00102 for (int i=0; i<numberOfMethods; ++i) 00103 { 00104 error = methods[i]->generateForecast(this, history, historycount, weight, debug); 00105 if (error<best_error) 00106 { 00107 best_error = error; 00108 best_method = i; 00109 } 00110 } 00111 } 00112 catch (...) 00113 { 00114 delete weight; 00115 throw; 00116 } 00117 delete weight; 00118 00119 // Apply the most appropriate forecasting method 00120 if (best_method >= 0) 00121 { 00122 if (debug) 00123 logger << getName() << ": chosen method: " << methods[best_method]->getName() << endl; 00124 methods[best_method]->applyForecast(this, buckets, bucketcount, debug); 00125 } 00126 } 00127 00128 00129 // 00130 // MOVING AVERAGE FORECAST 00131 // 00132 00133 00134 unsigned int Forecast::MovingAverage::defaultbuckets = 5; 00135 00136 00137 double Forecast::MovingAverage::generateForecast 00138 (Forecast* fcst, const double history[], unsigned int count, const double weight[], bool debug) 00139 { 00140 double error_smape = 0.0; 00141 00142 // Calculate the forecast and forecast error. 00143 for (unsigned int i = 1; i <= count; ++i) 00144 { 00145 double sum = 0.0; 00146 if (i > buckets) 00147 { 00148 for (unsigned int j = 0; j < buckets; ++j) 00149 sum += history[i-j-1]; 00150 avg = sum / buckets; 00151 } 00152 else 00153 { 00154 // For the first few values 00155 for (unsigned int j = 0; j < i; ++j) 00156 sum += history[i-j-1]; 00157 avg = sum / i; 00158 } 00159 if (i >= fcst->getForecastSkip() && i < count && avg + history[i] > ROUNDING_ERROR) 00160 error_smape += fabs(avg - history[i]) / (avg + history[i]) * weight[i]; 00161 } 00162 00163 // Echo the result 00164 if (debug) 00165 logger << (fcst ? fcst->getName() : "") << ": moving average : " 00166 << "smape " << error_smape 00167 << ", forecast " << avg << endl; 00168 return error_smape; 00169 } 00170 00171 00172 void Forecast::MovingAverage::applyForecast 00173 (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug) 00174 { 00175 // Loop over all buckets and set the forecast to a constant value 00176 if (avg < 0) return; 00177 for (unsigned int i = 1; i < bucketcount; ++i) 00178 forecast->setTotalQuantity( 00179 DateRange(buckets[i-1], buckets[i]), 00180 avg 00181 ); 00182 } 00183 00184 00185 // 00186 // SINGLE EXPONENTIAL FORECAST 00187 // 00188 00189 00190 double Forecast::SingleExponential::initial_alfa = 0.2; 00191 double Forecast::SingleExponential::min_alfa = 0.03; 00192 double Forecast::SingleExponential::max_alfa = 1.0; 00193 00194 00195 double Forecast::SingleExponential::generateForecast 00196 (Forecast* fcst, const double history[], unsigned int count, const double weight[], bool debug) 00197 { 00198 // Verify whether this is a valid forecast method. 00199 // - We need at least 5 buckets after the warmup period. 00200 if (count < fcst->getForecastSkip() + 5) 00201 return DBL_MAX; 00202 00203 unsigned int iteration = 1; 00204 bool upperboundarytested = false; 00205 bool lowerboundarytested = false; 00206 double error = 0.0, error_smape = 0.0, best_smape = 0.0, delta, df_dalfa_i, sum_11, sum_12; 00207 double best_error = DBL_MAX, best_alfa = initial_alfa, best_f_i = 0.0; 00208 for (; iteration <= Forecast::getForecastIterations(); ++iteration) 00209 { 00210 // Initialize variables 00211 df_dalfa_i = sum_11 = sum_12 = error_smape = error = 0.0; 00212 00213 // Initialize the iteration with the average of the first 3 values. 00214 f_i = (history[0] + history[1] + history[2]) / 3; 00215 00216 // Calculate the forecast and forecast error. 00217 // We also compute the sums required for the Marquardt optimization. 00218 for (unsigned long i = 1; i <= count; ++i) 00219 { 00220 df_dalfa_i = history[i-1] - f_i + (1 - alfa) * df_dalfa_i; 00221 f_i = history[i-1] * alfa + (1 - alfa) * f_i; 00222 if (i == count) break; 00223 sum_12 += df_dalfa_i * (history[i] - f_i) * weight[i]; 00224 sum_11 += df_dalfa_i * df_dalfa_i * weight[i]; 00225 if (i >= fcst->getForecastSkip()) 00226 { 00227 error += (f_i - history[i]) * (f_i - history[i]) * weight[i]; 00228 if (f_i + history[i] > ROUNDING_ERROR) 00229 error_smape += fabs(f_i - history[i]) / (f_i + history[i]) * weight[i]; 00230 } 00231 } 00232 00233 // Better than earlier iterations? 00234 if (error < best_error) 00235 { 00236 best_error = error; 00237 best_smape = error_smape; 00238 best_alfa = alfa; 00239 best_f_i = f_i; 00240 } 00241 00242 // Add Levenberg - Marquardt damping factor 00243 if (fabs(sum_11 + error / iteration) > ROUNDING_ERROR) 00244 sum_11 += error / iteration; 00245 00246 // Calculate a delta for the alfa parameter 00247 if (fabs(sum_11) < ROUNDING_ERROR) break; 00248 delta = sum_12 / sum_11; 00249 00250 // Stop when we are close enough and have tried hard enough 00251 if (fabs(delta) < ACCURACY && iteration > 3) break; 00252 00253 // New alfa 00254 alfa += delta; 00255 00256 // Stop when we repeatedly bounce against the limits. 00257 // Testing a limits once is allowed. 00258 if (alfa > max_alfa) 00259 { 00260 alfa = max_alfa; 00261 if (upperboundarytested) break; 00262 upperboundarytested = true; 00263 } 00264 else if (alfa < min_alfa) 00265 { 00266 alfa = min_alfa; 00267 if (lowerboundarytested) break; 00268 lowerboundarytested = true; 00269 } 00270 } 00271 00272 // Keep the best result 00273 f_i = best_f_i; 00274 00275 // Echo the result 00276 if (debug) 00277 logger << (fcst ? fcst->getName() : "") << ": single exponential : " 00278 << "alfa " << best_alfa 00279 << ", smape " << best_smape 00280 << ", " << iteration << " iterations" 00281 << ", forecast " << f_i << endl; 00282 return best_smape; 00283 } 00284 00285 00286 void Forecast::SingleExponential::applyForecast 00287 (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug) 00288 { 00289 // Loop over all buckets and set the forecast to a constant value 00290 if (f_i < 0) return; 00291 for (unsigned int i = 1; i < bucketcount; ++i) 00292 forecast->setTotalQuantity( 00293 DateRange(buckets[i-1], buckets[i]), 00294 f_i 00295 ); 00296 } 00297 00298 00299 // 00300 // DOUBLE EXPONENTIAL FORECAST 00301 // 00302 00303 double Forecast::DoubleExponential::initial_alfa = 0.2; 00304 double Forecast::DoubleExponential::min_alfa = 0.02; 00305 double Forecast::DoubleExponential::max_alfa = 1.0; 00306 double Forecast::DoubleExponential::initial_gamma = 0.2; 00307 double Forecast::DoubleExponential::min_gamma = 0.05; 00308 double Forecast::DoubleExponential::max_gamma = 1.0; 00309 double Forecast::DoubleExponential::dampenTrend = 0.8; 00310 00311 00312 double Forecast::DoubleExponential::generateForecast 00313 (Forecast* fcst, const double history[], unsigned int count, const double weight[], bool debug) 00314 { 00315 // Verify whether this is a valid forecast method. 00316 // - We need at least 5 buckets after the warmup period. 00317 if (count < fcst->getForecastSkip() + 5) 00318 return DBL_MAX; 00319 00320 // Define variables 00321 double error = 0.0, error_smape = 0.0, delta_alfa, delta_gamma, determinant; 00322 double constant_i_prev, trend_i_prev, d_constant_d_gamma_prev, 00323 d_constant_d_alfa_prev, d_constant_d_alfa, d_constant_d_gamma, 00324 d_trend_d_alfa, d_trend_d_gamma, d_forecast_d_alfa, d_forecast_d_gamma, 00325 sum11, sum12, sum22, sum13, sum23; 00326 double best_error = DBL_MAX, best_smape = 0, best_alfa = initial_alfa, 00327 best_gamma = initial_gamma, best_constant_i = 0.0, best_trend_i = 0.0; 00328 00329 // Iterations 00330 unsigned int iteration = 1, boundarytested = 0; 00331 for (; iteration <= Forecast::getForecastIterations(); ++iteration) 00332 { 00333 // Initialize variables 00334 error = error_smape = sum11 = sum12 = sum22 = sum13 = sum23 = 0.0; 00335 d_constant_d_alfa = d_constant_d_gamma = d_trend_d_alfa = d_trend_d_gamma = 0.0; 00336 d_forecast_d_alfa = d_forecast_d_gamma = 0.0; 00337 00338 // Initialize the iteration 00339 constant_i = (history[0] + history[1] + history[2]) / 3; 00340 trend_i = (history[3] - history[0]) / 3; 00341 00342 // Calculate the forecast and forecast error. 00343 // We also compute the sums required for the Marquardt optimization. 00344 for (unsigned long i = 1; i <= count; ++i) 00345 { 00346 constant_i_prev = constant_i; 00347 trend_i_prev = trend_i; 00348 constant_i = history[i-1] * alfa + (1 - alfa) * (constant_i_prev + trend_i_prev); 00349 trend_i = gamma * (constant_i - constant_i_prev) + (1 - gamma) * trend_i_prev; 00350 if (i == count) break; 00351 d_constant_d_gamma_prev = d_constant_d_gamma; 00352 d_constant_d_alfa_prev = d_constant_d_alfa; 00353 d_constant_d_alfa = history[i-1] - constant_i_prev - trend_i_prev 00354 + (1 - alfa) * d_forecast_d_alfa; 00355 d_constant_d_gamma = (1 - alfa) * d_forecast_d_gamma; 00356 d_trend_d_alfa = gamma * (d_constant_d_alfa - d_constant_d_alfa_prev) 00357 + (1 - gamma) * d_trend_d_alfa; 00358 d_trend_d_gamma = constant_i - constant_i_prev - trend_i_prev 00359 + gamma * (d_constant_d_gamma - d_constant_d_gamma_prev) 00360 + (1 - gamma) * d_trend_d_gamma; 00361 d_forecast_d_alfa = d_constant_d_alfa + d_trend_d_alfa; 00362 d_forecast_d_gamma = d_constant_d_gamma + d_trend_d_gamma; 00363 sum11 += weight[i] * d_forecast_d_alfa * d_forecast_d_alfa; 00364 sum12 += weight[i] * d_forecast_d_alfa * d_forecast_d_gamma; 00365 sum22 += weight[i] * d_forecast_d_gamma * d_forecast_d_gamma; 00366 sum13 += weight[i] * d_forecast_d_alfa * (history[i] - constant_i - trend_i); 00367 sum23 += weight[i] * d_forecast_d_gamma * (history[i] - constant_i - trend_i); 00368 if (i >= fcst->getForecastSkip()) // Don't measure during the warmup period 00369 { 00370 error += (constant_i + trend_i - history[i]) * (constant_i + trend_i - history[i]) * weight[i]; 00371 if (constant_i + trend_i + history[i] > ROUNDING_ERROR) 00372 error_smape += fabs(constant_i + trend_i - history[i]) / (constant_i + trend_i + history[i]) * weight[i]; 00373 } 00374 } 00375 00376 // Better than earlier iterations? 00377 if (error < best_error) 00378 { 00379 best_error = error; 00380 best_smape = error_smape; 00381 best_alfa = alfa; 00382 best_gamma = gamma; 00383 best_constant_i = constant_i; 00384 best_trend_i = trend_i; 00385 } 00386 00387 // Add Levenberg - Marquardt damping factor 00388 sum11 += error / iteration; 00389 sum22 += error / iteration; 00390 00391 // Calculate a delta for the alfa and gamma parameters 00392 determinant = sum11 * sum22 - sum12 * sum12; 00393 if (fabs(determinant) < ROUNDING_ERROR) 00394 { 00395 // Almost singular matrix. Try without the damping factor. 00396 sum11 -= error / iteration; 00397 sum22 -= error / iteration; 00398 determinant = sum11 * sum22 - sum12 * sum12; 00399 if (fabs(determinant) < ROUNDING_ERROR) 00400 // Still singular - stop iterations here 00401 break; 00402 } 00403 delta_alfa = (sum13 * sum22 - sum23 * sum12) / determinant; 00404 delta_gamma = (sum23 * sum11 - sum13 * sum12) / determinant; 00405 00406 // Stop when we are close enough and have tried hard enough 00407 if (fabs(delta_alfa) + fabs(delta_gamma) < 2 * ACCURACY && iteration > 3) 00408 break; 00409 00410 // New values for the next iteration 00411 alfa += delta_alfa; 00412 gamma += delta_gamma; 00413 00414 // Limit the parameters in their allowed range. 00415 if (alfa > max_alfa) 00416 alfa = max_alfa; 00417 else if (alfa < min_alfa) 00418 alfa = min_alfa; 00419 if (gamma > max_gamma) 00420 gamma = max_gamma; 00421 else if (gamma < min_gamma) 00422 gamma = min_gamma; 00423 00424 // Verify repeated running with both parameters at the boundary 00425 if ((gamma == min_gamma || gamma == max_gamma) 00426 && (alfa == min_alfa || alfa == max_alfa)) 00427 { 00428 if (boundarytested++ > 5) break; 00429 } 00430 } 00431 00432 // Keep the best result 00433 constant_i = best_constant_i; 00434 trend_i = best_trend_i; 00435 00436 // Echo the result 00437 if (debug) 00438 logger << (fcst ? fcst->getName() : "") << ": double exponential : " 00439 << "alfa " << best_alfa 00440 << ", gamma " << best_gamma 00441 << ", smape " << best_smape 00442 << ", " << iteration << " iterations" 00443 << ", constant " << constant_i 00444 << ", trend " << trend_i 00445 << ", forecast " << (trend_i + constant_i) << endl; 00446 return best_smape; 00447 } 00448 00449 00450 void Forecast::DoubleExponential::applyForecast 00451 (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug) 00452 { 00453 // Loop over all buckets and set the forecast to a linearly changing value 00454 for (unsigned int i = 1; i < bucketcount; ++i) 00455 { 00456 constant_i += trend_i; 00457 trend_i *= dampenTrend; // Reduce slope in the future 00458 if (constant_i > 0) 00459 forecast->setTotalQuantity( 00460 DateRange(buckets[i-1], buckets[i]), 00461 constant_i 00462 ); 00463 } 00464 } 00465 00466 00467 // 00468 // SEASONAL FORECAST 00469 // 00470 00471 unsigned int Forecast::Seasonal::min_period = 2; 00472 unsigned int Forecast::Seasonal::max_period = 14; 00473 double Forecast::Seasonal::dampenTrend = 0.8; 00474 double Forecast::Seasonal::initial_alfa = 0.2; 00475 double Forecast::Seasonal::min_alfa = 0.02; 00476 double Forecast::Seasonal::max_alfa = 1.0; 00477 double Forecast::Seasonal::initial_beta = 0.2; 00478 double Forecast::Seasonal::min_beta = 0.2; 00479 double Forecast::Seasonal::max_beta = 1.0; 00480 double Forecast::Seasonal::initial_gamma = 0.3; 00481 double Forecast::Seasonal::min_gamma = 0.05; 00482 double Forecast::Seasonal::max_gamma = 1.0; 00483 00484 00485 void Forecast::Seasonal::detectCycle(const double history[], unsigned int count) 00486 { 00487 // We need at least 2 cycles 00488 if (count < min_period*2) return; 00489 00490 // Compute the average value 00491 double average = 0.0; 00492 for (unsigned int i = 0; i < count; ++i) 00493 average += history[i]; 00494 average /= count; 00495 00496 // Compute variance 00497 double variance = 0.0; 00498 for (unsigned int i = 0; i < count; ++i) 00499 variance += (history[i]-average)*(history[i]-average); 00500 variance /= count; 00501 00502 // Compute autocorrelation for different periods 00503 double prev = 10.0; 00504 double prevprev = 10.0; 00505 for (unsigned short p = min_period; p <= max_period && p < count/2; ++p) 00506 { 00507 // Compute correlation 00508 double correlation = 0.0; 00509 for (unsigned int i = p; i < count; ++i) 00510 correlation += (history[i]-average)*(history[i-p]-average); 00511 correlation /= count - p; 00512 correlation /= variance; 00513 // Detect cycles if we find a period with a big autocorrelation which 00514 // is significantly larger than the adjacent periods. 00515 if ((p > min_period + 1 && prev > prevprev*1.1) && prev > correlation*1.1 && prev > 0.3) 00516 { 00517 period = p - 1; 00518 return; 00519 } 00520 prevprev = prev; 00521 prev = correlation; 00522 } 00523 } 00524 00525 00526 double Forecast::Seasonal::generateForecast 00527 (Forecast* fcst, const double history[], unsigned int count, const double weight[], bool debug) 00528 { 00529 // Check for seasonal cycles 00530 detectCycle(history, count); 00531 00532 // Return if no seasonality is found 00533 if (!period) return DBL_MAX; 00534 00535 // Define variables 00536 double error = 0.0, error_smape = 0.0, delta_alfa, delta_beta, delta_gamma; 00537 double sum11, sum12, sum13, sum14, sum22, sum23, sum24, sum33, sum34; 00538 double best_error = DBL_MAX, best_smape = 0, best_alfa = initial_alfa, 00539 best_beta = initial_beta, best_gamma = initial_gamma; 00540 S_i = new double[period]; 00541 00542 // Iterations 00543 double L_i_prev; 00544 unsigned int iteration = 1, boundarytested = 0; 00545 for (; iteration <= Forecast::getForecastIterations(); ++iteration) 00546 { 00547 // Initialize variables 00548 error = error_smape = sum11 = sum12 = sum13 = sum14 = 0.0; 00549 sum22 = sum23 = sum24 = sum33 = sum34 = 0.0; 00550 00551 // Initialize the iteration 00552 // L_i = average over first cycle 00553 // T_i = average delta measured in second cycle 00554 // S_i[index] = actual divided by average in first cycle 00555 L_i = 0.0; 00556 for (unsigned long i = 0; i < period; ++i) 00557 L_i += history[i]; 00558 L_i /= period; 00559 T_i = 0.0; 00560 for (unsigned long i = 0; i < period; ++i) 00561 { 00562 T_i += history[i+period] - history[i]; 00563 S_i[i] = history[i] / L_i; 00564 } 00565 T_i /= period * period; 00566 00567 // Calculate the forecast and forecast error. 00568 // We also compute the sums required for the Marquardt optimization. 00569 cycleindex = 0; 00570 for (unsigned long i = period; i <= count; ++i) 00571 { 00572 L_i_prev = L_i; 00573 if (S_i[cycleindex] > ROUNDING_ERROR) 00574 L_i = alfa * history[i-1] / S_i[cycleindex] + (1 - alfa) * (L_i + T_i); 00575 else 00576 L_i = (1 - alfa) * (L_i + T_i); 00577 T_i = beta * (L_i - L_i_prev) + (1 - beta) * T_i; 00578 S_i[cycleindex] = gamma * history[i-1] / L_i + (1 - gamma) * S_i[cycleindex]; 00579 if (i == count) break; 00580 if (i >= fcst->getForecastSkip()) // Don't measure during the warmup period 00581 { 00582 double fcst = (L_i + T_i) * S_i[cycleindex]; 00583 error += (fcst - history[i]) * (fcst - history[i]) * weight[i]; 00584 if (fcst + history[i] > ROUNDING_ERROR) 00585 error_smape += fabs(fcst - history[i]) / (fcst + history[i]) * weight[i]; 00586 } 00587 if (++cycleindex >= period) cycleindex = 0; 00588 } 00589 00590 // Better than earlier iterations? 00591 best_smape = error_smape; 00592 if (error < best_error) 00593 { 00594 best_error = error; 00595 best_smape = error_smape; 00596 best_alfa = alfa; 00597 best_beta = beta; 00598 best_gamma = gamma; 00599 } 00600 break; // @todo no iterations yet to tune the seasonal parameters 00601 00602 // Add Levenberg - Marquardt damping factor 00603 sum11 += error / iteration; 00604 sum22 += error / iteration; 00605 sum33 += error / iteration; 00606 00607 // Calculate a delta for the alfa, beta and gamma parameters. 00608 // We're using Cramer's rule to solve a set of linear equations. 00609 double det = determinant(sum11, sum12, sum13, 00610 sum12, sum22, sum23, 00611 sum13, sum23, sum33); 00612 if (fabs(det) < ROUNDING_ERROR) 00613 { 00614 // Almost singular matrix. Try without the damping factor. 00615 sum11 -= error / iteration; 00616 sum22 -= error / iteration; 00617 sum33 -= error / iteration; 00618 det = determinant(sum11, sum12, sum13, 00619 sum12, sum22, sum23, 00620 sum13, sum23, sum33); 00621 if (fabs(det) < ROUNDING_ERROR) 00622 // Still singular - stop iterations here 00623 break; 00624 } 00625 delta_alfa = determinant(sum14, sum12, sum13, 00626 sum24, sum22, sum23, 00627 sum34, sum23, sum33) / det; 00628 delta_beta = determinant(sum11, sum14, sum13, 00629 sum12, sum24, sum23, 00630 sum13, sum34, sum33) / det; 00631 delta_gamma = determinant(sum11, sum12, sum14, 00632 sum12, sum22, sum24, 00633 sum13, sum23, sum34) / det; 00634 00635 // Stop when we are close enough and have tried hard enough 00636 if (fabs(delta_alfa) + fabs(delta_beta) + fabs(delta_gamma) < 3 * ACCURACY 00637 && iteration > 3) 00638 break; 00639 00640 // New values for the next iteration 00641 alfa += delta_alfa; 00642 alfa += delta_beta; 00643 gamma += delta_gamma; 00644 00645 // Limit the parameters in their allowed range. 00646 if (alfa > max_alfa) 00647 alfa = max_alfa; 00648 else if (alfa < min_alfa) 00649 alfa = min_alfa; 00650 if (beta > max_beta) 00651 beta = max_beta; 00652 else if (beta < min_beta) 00653 beta = min_beta; 00654 if (gamma > max_gamma) 00655 gamma = max_gamma; 00656 else if (gamma < min_gamma) 00657 gamma = min_gamma; 00658 00659 // Verify repeated running with any parameters at the boundary 00660 if (gamma == min_gamma || gamma == max_gamma || 00661 beta == min_beta || beta == max_beta || 00662 alfa == min_alfa || alfa == max_alfa) 00663 { 00664 if (boundarytested++ > 5) break; 00665 } 00666 } 00667 00668 if (period > fcst->getForecastSkip()) 00669 // Correction on the error: we counted less buckets. We now 00670 // proportionally increase the error to account for this and have a 00671 // value that can be compared with the other forecast methods. 00672 best_smape *= (count - fcst->getForecastSkip()) / (count - period); 00673 00674 // Keep the best result 00675 00676 // Echo the result 00677 if (debug) 00678 logger << (fcst ? fcst->getName() : "") << ": seasonal : " 00679 << "alfa " << best_alfa 00680 << ", beta " << best_beta 00681 << ", gamma " << best_gamma 00682 << ", smape " << best_smape 00683 << ", " << iteration << " iterations" 00684 << ", period " << period 00685 << ", constant " << L_i 00686 << ", trend " << T_i 00687 << ", forecast " << (L_i + T_i) * S_i[count % period] 00688 << endl; 00689 return best_smape; 00690 } 00691 00692 00693 void Forecast::Seasonal::applyForecast 00694 (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug) 00695 { 00696 // Loop over all buckets and set the forecast to a linearly changing value 00697 for (unsigned int i = 1; i < bucketcount; ++i) 00698 { 00699 L_i += T_i; 00700 T_i *= dampenTrend; // Reduce slope in the future 00701 if (L_i * S_i[cycleindex] > 0) 00702 forecast->setTotalQuantity( 00703 DateRange(buckets[i-1], buckets[i]), 00704 L_i * S_i[cycleindex] 00705 ); 00706 if (++cycleindex >= period) cycleindex = 0; 00707 } 00708 } 00709 00710 00711 // 00712 // CROSTON'S FORECAST METHOD 00713 // 00714 00715 00716 double Forecast::Croston::initial_alfa = 0.1; 00717 double Forecast::Croston::min_alfa = 0.03; 00718 double Forecast::Croston::max_alfa = 1.0; 00719 double Forecast::Croston::min_intermittence = 0.33; 00720 00721 00722 double Forecast::Croston::generateForecast 00723 (Forecast* fcst, const double history[], unsigned int count, const double weight[], bool debug) 00724 { 00725 unsigned int iteration = 1; 00726 bool upperboundarytested = false; 00727 bool lowerboundarytested = false; 00728 double error = 0.0, error_smape = 0.0, best_smape = 0.0, delta; 00729 double q_i, p_i, d_p_i, d_q_i, d_f_i, sum1, sum2; 00730 double best_error = DBL_MAX, best_alfa = initial_alfa, best_f_i = 0.0; 00731 unsigned int between_demands = 1; 00732 for (; iteration <= Forecast::getForecastIterations(); ++iteration) 00733 { 00734 // Initialize variables 00735 error_smape = error = d_p_i = d_q_i = d_f_i = sum1 = sum2 = 0.0; 00736 00737 // Initialize the iteration. 00738 q_i = f_i = history[0]; 00739 p_i = 0; 00740 00741 // Calculate the forecast and forecast error. 00742 // We also compute the sums required for the Marquardt optimization. 00743 for (unsigned long i = 1; i <= count; ++i) 00744 { 00745 if (history[i-1]) 00746 { 00747 // Non-zero bucket 00748 d_p_i = between_demands - p_i + (1 - alfa) * d_p_i; 00749 d_q_i = history[i-1] - q_i + (1 - alfa) * d_q_i; 00750 q_i = alfa * history[i-1] + (1 - alfa) * q_i; 00751 p_i = alfa * between_demands + (1 - alfa) * q_i; 00752 f_i = q_i / p_i; 00753 d_f_i = (d_q_i - d_p_i * q_i / p_i) / p_i; 00754 between_demands = 1; 00755 } 00756 else 00757 ++between_demands; 00758 if (i == count) break; 00759 sum1 += weight[i] * d_f_i * (history[i] - f_i); 00760 sum2 += weight[i] * d_f_i * d_f_i; 00761 if (i >= fcst->getForecastSkip() && p_i > 0) 00762 { 00763 error += (f_i - history[i]) * (f_i - history[i]) * weight[i]; 00764 if (f_i + history[i] > ROUNDING_ERROR) 00765 error_smape += fabs(f_i - history[i]) / (f_i + history[i]) * weight[i]; 00766 } 00767 } 00768 00769 // Better than earlier iterations? 00770 if (error < best_error) 00771 { 00772 best_error = error; 00773 best_smape = error_smape; 00774 best_alfa = alfa; 00775 best_f_i = f_i; 00776 } 00777 00778 // Add Levenberg - Marquardt damping factor 00779 if (fabs(sum2 + error / iteration) > ROUNDING_ERROR) 00780 sum2 += error / iteration; 00781 00782 // Calculate a delta for the alfa parameter 00783 if (fabs(sum2) < ROUNDING_ERROR) break; 00784 delta = sum1 / sum2; 00785 00786 // Stop when we are close enough and have tried hard enough 00787 if (fabs(delta) < ACCURACY && iteration > 3) break; 00788 00789 // New alfa 00790 alfa += delta; 00791 00792 // Stop when we repeatedly bounce against the limits. 00793 // Testing a limits once is allowed. 00794 if (alfa > max_alfa) 00795 { 00796 alfa = max_alfa; 00797 if (upperboundarytested) break; 00798 upperboundarytested = true; 00799 } 00800 else if (alfa < min_alfa) 00801 { 00802 alfa = min_alfa; 00803 if (lowerboundarytested) break; 00804 lowerboundarytested = true; 00805 } 00806 } 00807 00808 // Keep the best result 00809 f_i = best_f_i; 00810 alfa = best_alfa; 00811 00812 // Echo the result 00813 if (debug) 00814 logger << (fcst ? fcst->getName() : "") << ": croston : " 00815 << "alfa " << best_alfa 00816 << ", smape " << best_smape 00817 << ", " << iteration << " iterations" 00818 << ", forecast " << f_i << endl; 00819 return best_smape; 00820 } 00821 00822 00823 void Forecast::Croston::applyForecast 00824 (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug) 00825 { 00826 // Loop over all buckets and set the forecast to a constant value 00827 if (f_i < 0) return; 00828 for (unsigned int i = 1; i < bucketcount; ++i) 00829 forecast->setTotalQuantity( 00830 DateRange(buckets[i-1], buckets[i]), 00831 f_i 00832 ); 00833 } 00834 00835 } // end namespace