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