Holt-Winters calculation improvement

2 files changed, 123 insertions, 10 deletions

diff --git a/example/ndpiReader.c b/example/ndpiReader.c
index e268b0198..207f6a2ad 100644
--- a/example/ndpiReader.c
+++ b/example/ndpiReader.c
@@ -3787,6 +3787,114 @@ void hwUnitTest() {
 
 /* *********************************************** */
 
+void hwUnitTest2() {
+  struct ndpi_hw_struct hw;
+  u_int num_learning_points = 12;
+  u_int8_t trace = 1;
+  double v[] = {
+    31.908466339111,
+    87.339714050293,
+    173.47660827637,
+    213.92568969727,
+    223.32124328613,
+    230.60134887695,
+    238.09457397461,
+    245.8137512207,
+    251.09228515625,
+    251.09228515625,
+    259.21997070312,
+    261.98754882812,
+    264.78540039062,
+    264.78540039062,
+    270.47451782227,
+    273.3671875,
+    288.34222412109,
+    288.34222412109,
+    304.24795532227,
+    304.24795532227,
+    350.92227172852,
+    384.54431152344,
+    423.25942993164,
+    439.43322753906,
+    445.05981445312,
+    445.05981445312,
+    445.05981445312,
+    445.05981445312
+  };
+  u_int i, num = sizeof(v) / sizeof(double);
+  float alpha = 0.5, beta = 0.5, gamma = 0.1;
+
+  assert(ndpi_hw_init(&hw, num_learning_points, 0 /* 0=multiplicative, 1=additive */, alpha, beta, gamma, 0.05) == 0);
+  
+  if(trace)
+    printf("\nHolt-Winters [alpha: %.1f][beta: %.1f][gamma: %.1f]\n", alpha, beta, gamma);
+  
+  for(i=0; i<num; i++) {
+    double prediction, confidence_band;
+    double lower, upper;
+    int rc = ndpi_hw_add_value(&hw, v[i], &prediction, &confidence_band);
+    
+    lower = prediction - confidence_band, upper = prediction + confidence_band;
+    
+    if(trace)
+      printf("%2u)\t%12.3f\t%.3f\t%12.3f\t%12.3f\t %s [%.3f]\n", i, v[i], prediction, lower, upper,
+	     ((rc == 0) || ((v[i] >= lower) && (v[i] <= upper))) ? "OK" : "ANOMALY",
+	     confidence_band);
+  }
+  
+  ndpi_hw_free(&hw);
+}
+
+/* *********************************************** */
+
+void hwUnitTest3() {
+  struct ndpi_hw_struct hw;
+  u_int num_learning_points = 3;
+  u_int8_t trace = 1;
+  double v[] = {
+    10,
+    14,
+    8,
+    25,
+    16,
+    22,
+    14,
+    35,
+    15,
+    27,
+    18,
+    40,
+    28,
+    40,
+    25,
+    65,
+  };
+  u_int i, num = sizeof(v) / sizeof(double);
+  float alpha = 0.5, beta = 0.5, gamma = 0.1;
+  assert(ndpi_hw_init(&hw, num_learning_points, 0 /* 0=multiplicative, 1=additive */, alpha, beta, gamma, 0.05) == 0);
+  
+  if(trace)
+    printf("\nHolt-Winters [alpha: %.1f][beta: %.1f][gamma: %.1f]\n", alpha, beta, gamma);
+  
+  for(i=0; i<num; i++) {
+    double prediction, confidence_band;
+    double lower, upper;
+    int rc = ndpi_hw_add_value(&hw, v[i], &prediction, &confidence_band);
+    
+    lower = prediction - confidence_band, upper = prediction + confidence_band;
+    
+    if(trace)
+      printf("%2u)\t%12.3f\t%.3f\t%12.3f\t%12.3f\t %s [%.3f]\n",
+	     i, v[i], prediction, lower, upper,
+	     ((rc == 0) || ((v[i] >= lower) && (v[i] <= upper))) ? "OK" : "ANOMALY",
+	     confidence_band);
+  }
+  
+  ndpi_hw_free(&hw);
+}
+
+/* *********************************************** */
+
 void jitterUnitTest() {
   struct ndpi_jitter_struct jitter;
   float v[] = { 10, 14, 8, 25, 16, 22, 14, 35, 15, 27, 218, 40, 28, 40, 25, 65 };
@@ -3827,7 +3935,9 @@ int orginal_main(int argc, char **argv) {
     ndpi_info_mod = ndpi_init_detection_module(ndpi_no_prefs);
 
     if(ndpi_info_mod == NULL) return -1;
-
+    
+    // hwUnitTest2();
+    
     /* Internal checks */    
     // binUnitTest();
     hwUnitTest();
diff --git a/src/lib/ndpi_analyze.c b/src/lib/ndpi_analyze.c
index e142f914c..c46804ead 100644
--- a/src/lib/ndpi_analyze.c
+++ b/src/lib/ndpi_analyze.c
@@ -969,7 +969,7 @@ int ndpi_hw_init(struct ndpi_hw_struct *hw,
 		 double alpha, double beta, double gamma, float significance) {
   memset(hw, 0, sizeof(struct ndpi_hw_struct));
 
-  hw->params.num_season_periods = num_periods;
+  hw->params.num_season_periods = num_periods + 1;
   hw->params.alpha      = alpha;
   hw->params.beta       = beta;
   hw->params.gamma      = gamma;
@@ -978,6 +978,8 @@ int ndpi_hw_init(struct ndpi_hw_struct *hw,
   if((significance < 0) || (significance > 1)) significance = 0.05;
   hw->params.ro         = ndpi_normal_cdf_inverse(1 - (significance / 2.));
 
+  hw->params.ro = 1.95996398454005;
+  
   if((hw->y = (u_int32_t*)ndpi_calloc(hw->params.num_season_periods, sizeof(u_int32_t))) == NULL)
     return(-1);
 
@@ -1038,14 +1040,12 @@ int ndpi_hw_add_value(struct ndpi_hw_struct *hw, const u_int32_t _value, double
     double prev_u = hw->u;
     double prev_v = hw->v;
     double prev_s = hw->s[idx];
-    double value  = (double)_value;;
-    double error;
-
-    hw->u      = ((hw->params.alpha * value) / prev_s) + (1 - hw->params.alpha) * (hw->u +  hw->v);
-    hw->v      = (hw->params.beta * (hw->u - prev_u))  + ((1 - hw->params.beta) * hw->v);
-    hw->s[idx] = (hw->params.gamma * (value / hw->u))  + ((1 - hw->params.gamma) * prev_s);
+    double value  = (double)_value;
+    double sq, error;
 
-    hw->num_values++, idx = (idx + 1) % hw->params.num_season_periods;
+    hw->u      = ((hw->params.alpha * value) / prev_s)  + ( 1 - hw->params.alpha) * (hw->u + hw->v);
+    hw->v      = (hw->params.beta   * (hw->u - prev_u)) + ((1 - hw->params.beta ) * hw->v);
+    hw->s[idx] = (hw->params.gamma  * (value / hw->u))  + ((1 - hw->params.gamma) * prev_s);
 
     if(hw->params.use_hw_additive_seasonal)
       *forecast = (prev_u + prev_v) + prev_s;
@@ -1054,8 +1054,11 @@ int ndpi_hw_add_value(struct ndpi_hw_struct *hw, const u_int32_t _value, double
 
     error                 = value - *forecast;
     hw->sum_square_error += error * error;
-    *confidence_band      = hw->params.ro * (sqrt(hw->sum_square_error / (hw->num_values - hw->params.num_season_periods)));
+    sq = sqrt(hw->sum_square_error / (hw->num_values + 2 - hw->params.num_season_periods));
+    *confidence_band      = hw->params.ro * sq;
 
+    hw->num_values++, idx = (idx + 1) % hw->params.num_season_periods;
+    
     return(1); /* We're in business: forecast is meaningful now */
   }
 }