From 13c63f495ce7c0a47ceef6ce7af4559cd1aaa6f4 Mon Sep 17 00:00:00 2001
From: Ruben Dashyan <ruben.dashyan@numworks.com>
Date: Tue, 28 May 2019 12:17:47 +0200
Subject: [PATCH] [apps/regression/model/exponential_model] Explicit fit

---
 apps/regression/model/exponential_model.cpp | 39 +++++++++++++++++++++
 apps/regression/model/exponential_model.h   |  1 +
 2 files changed, 40 insertions(+)

diff --git a/apps/regression/model/exponential_model.cpp b/apps/regression/model/exponential_model.cpp
index e29bdb047..ee7ea3c60 100644
--- a/apps/regression/model/exponential_model.cpp
+++ b/apps/regression/model/exponential_model.cpp
@@ -1,4 +1,5 @@
 #include "exponential_model.h"
+#include "../store.h"
 #include <math.h>
 #include <assert.h>
 #include <poincare/code_point_layout.h>
@@ -46,6 +47,44 @@ double ExponentialModel::levelSet(double * modelCoefficients, double xMin, doubl
   return log(y/a)/b;
 }
 
+void ExponentialModel::fit(Store * store, int series, double * modelCoefficients, Poincare::Context * context) {
+  /* By the change of variable z=ln(y), the equation y=a*exp(b*x) becomes
+   * z=c*x+d with c=b and d=ln(a). Although that change of variable does not
+   * preserve the regression error function, it turns an exponential regression
+   * problem into a linear one and we consider that the solution of the latter
+   * is good enough for our purpose.
+   * That being said, one should check that the y values are all positive. (If
+   * the y values are all negative, one may replace each of them by its
+   * opposite. In the case where y values happen to be zero or of opposite
+   * sign, we call the base class method as a fallback. */
+  double sumOfX = 0;
+  double sumOfY = 0;
+  double sumOfXX = 0;
+  double sumOfXY = 0;
+  const int numberOfPoints = store->numberOfPairsOfSeries(series);
+  const int sign = store->get(series, 1, 0) > 0 ? 1 : -1;
+  for (int p = 0; p < numberOfPoints; p++) {
+    const double x = store->get(series, 0, p);
+    const double z = store->get(series, 1, p) * sign;
+    if (z <= 0) {
+      return Model::fit(store, series, modelCoefficients, context);
+    }
+    const double y = log(z);
+    sumOfX += x;
+    sumOfY += y;
+    sumOfXX += x*x;
+    sumOfXY += x*y;
+  }
+  const double meanOfX = sumOfX / numberOfPoints;
+  const double meanOfY = sumOfY / numberOfPoints;
+  const double meanOfXX = sumOfXX / numberOfPoints;
+  const double meanOfXY = sumOfXY / numberOfPoints;
+  const double variance = meanOfXX - meanOfX * meanOfX;
+  const double covariance = meanOfXY - meanOfX * meanOfY;
+  modelCoefficients[1] = covariance / variance;
+  modelCoefficients[0] = sign * exp(meanOfY - modelCoefficients[1] * meanOfX);
+}
+
 double ExponentialModel::partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const {
   double a = modelCoefficients[0];
   double b = modelCoefficients[1];
diff --git a/apps/regression/model/exponential_model.h b/apps/regression/model/exponential_model.h
index 7d373b473..d7293ce35 100644
--- a/apps/regression/model/exponential_model.h
+++ b/apps/regression/model/exponential_model.h
@@ -12,6 +12,7 @@ public:
   I18n::Message formulaMessage() const override { return I18n::Message::ExponentialRegressionFormula; }
   double evaluate(double * modelCoefficients, double x) const override;
   double levelSet(double * modelCoefficients, double xMin, double step, double xMax, double y, Poincare::Context * context) override;
+  void fit(Store * store, int series, double * modelCoefficients, Poincare::Context * context) override;
   double partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const override;
   int numberOfCoefficients() const override { return 2; }
   int bannerLinesCount() const override { return 2; }