Merge commit '39f0ad12e577d8a7ad0c95f8dcfd1d01316185aa'

ARPA-SIMC · Jan 19, 2024 · 5df9a12 · 5df9a12
2 parents 74427a4 + 39f0ad1
commit 5df9a12
Show file tree

Hide file tree

Showing 4 changed files with 107 additions and 0 deletions.
diff --git a/agrolib/mathFunctions/furtherMathFunctions.cpp b/agrolib/mathFunctions/furtherMathFunctions.cpp
@@ -1000,6 +1000,36 @@ namespace interpolation
         return weighted_r_squared;
     }
 
+    double computeWeighted_StandardError(const std::vector<double>& observed, const std::vector<double>& predicted, const std::vector<double>& weights, int nrPredictors)
+    {
+        // This function computes the weighted R-squared (coefficient of determination)
+        double sum_weighted_squared_residuals = 0.0;
+        double sum_weighted_squared_total = 0.0;
+        double weighted_mean_observed = 0.0;
+
+        // Calculate the weighted mean of the observed values
+        double sum_weights = 0.0;
+        for (int i = 0; i < observed.size(); i++)
+        {
+            weighted_mean_observed += observed[i] * weights[i];
+            sum_weights += weights[i];
+        }
+        weighted_mean_observed /= sum_weights;
+
+        // Calculate the sums needed for weighted R-squared calculation
+        for (int i = 0; i < observed.size(); i++)
+        {
+            double weighted_residual = weights[i] * (observed[i] - predicted[i]);
+            sum_weighted_squared_residuals += weighted_residual * weighted_residual;
+        }
+        double standardError;
+        if (observed.size() > (nrPredictors+1))
+            standardError = sqrt(sum_weighted_squared_residuals/(observed.size()-nrPredictors-1));
+        else
+            standardError = sqrt(sum_weighted_squared_residuals/(observed.size()-1));
+        return 0;
+    }
+
 
     int bestFittingMarquardt_nDimension(double (*func)(std::vector<std::function<double(double, std::vector<double>&)>>&, std::vector<double>& , std::vector <std::vector <double>>&),
                                         std::vector<std::function<double(double, std::vector<double>&)>>& myFunc,

diff --git a/agrolib/mathFunctions/furtherMathFunctions.h b/agrolib/mathFunctions/furtherMathFunctions.h
@@ -91,6 +91,7 @@ enum estimatedFunction {FUNCTION_CODE_SPHERICAL, FUNCTION_CODE_LINEAR, FUNCTION_
 
         double computeR2(const std::vector<double>& obs, const std::vector<double>& sim);
         double computeWeighted_R2(const std::vector<double>& observed, const std::vector<double>& predicted, const std::vector<double>& weights);
+        double computeWeighted_StandardError(const std::vector<double>& observed, const std::vector<double>& predicted, const std::vector<double>& weights, int nrPredictors);
         double weightedVariance(const std::vector<double>& data, const std::vector<double>& weights);
         int bestFittingMarquardt_nDimension(double (*func)(std::vector<std::function<double (double, std::vector<double> &)> > &, std::vector<double> &, std::vector<std::vector<double>> &),
                                         std::vector<std::function<double (double, std::vector<double> &)> >& myFunc,

diff --git a/agrolib/mathFunctions/statistics.cpp b/agrolib/mathFunctions/statistics.cpp
@@ -499,6 +499,81 @@ namespace statistics
         free(roots);
     }
 
+    void weightedMultiRegressionLinearWithoutConstantTerm(const std::vector <std::vector <float>> &x, std::vector <float> &y, const std::vector <float> &weight, long nrItems,std::vector <float> &m, int nrPredictors)
+    {
+        double** XT = (double**)calloc(nrPredictors, sizeof(double*));
+        double** X = (double**)calloc(nrItems, sizeof(double*));
+        double** X2 = (double**)calloc(nrPredictors, sizeof(double*));
+        double** X2Inverse = (double**)calloc(nrPredictors, sizeof(double*));
+        for (int i=0;i<nrPredictors;i++)
+        {
+            XT[i]= (double*)calloc(nrItems, sizeof(double));
+            X2[i]= (double*)calloc(nrItems, sizeof(double));
+            X2Inverse[i]= (double*)calloc(nrItems, sizeof(double));
+        }
+        for (int i=0;i<nrItems;i++)
+        {
+            X[i]= (double*)calloc(nrPredictors, sizeof(double));
+        }
+        for (int j=0;j<nrPredictors;j++)
+        {
+            for (int i=0;i<nrItems;i++)
+            {
+                    XT[j][i] = x[i][j];
+            }
+        }
+        matricial::transposedMatrix(XT,nrPredictors,nrItems,X);
+        for (int j=0;j<nrPredictors;j++)
+        {
+            for (int i =0; i<nrItems; i++)
+            {
+                X[i][j]= 1./weight[i]*X[i][j];
+            }
+        }
+        matricial::matrixProduct(XT,X,nrItems,nrPredictors,nrPredictors,nrItems,X2);
+        matricial::inverse(X2,X2Inverse,nrPredictors);
+        //matricial::matrixProduct(X2Inverse,XT,nrPredictors+1,nrPredictors+1,nrItems,nrPredictors+1,X);
+        for (int i=0;i<nrItems;i++)
+        {
+            y[i] /= weight[i];
+        }
+        double* roots = (double*)calloc(nrPredictors, sizeof(double));
+        for (int j=0;j<nrPredictors;j++)
+        {
+            roots[j]=0;
+            for (int i=0;i<nrItems;i++)
+            {
+                roots[j] += (XT[j][i]*y[i]);
+            }
+        }
+        for (int j=0;j<nrPredictors;j++)
+        {
+            m[j]=0;
+        }
+        for (int j=0;j<nrPredictors;j++)
+        {
+            for (int i=0;i<nrPredictors;i++)
+            {
+                m[j] += float(X2Inverse[j][i]*roots[i]);
+            }
+        }
+        for (int j=0;j<nrPredictors;j++)
+        {
+            free(XT[j]);
+            free(X2[j]);
+            free(X2Inverse[j]);
+        }
+        for (int i=0;i<nrItems;i++)
+        {
+            free(X[i]);
+        }
+        free(X);
+        free(XT);
+        free(X2);
+        free(X2Inverse);
+        free(roots);
+    }
+
     void multiRegressionLinear(float** x,  float* y, long nrItems,float* q,float* m, int nrPredictors)
     {
         double* SUMx;

diff --git a/agrolib/mathFunctions/statistics.h b/agrolib/mathFunctions/statistics.h
@@ -33,6 +33,7 @@
         float linearInterpolation(float x1, float y1, float x2, float y2, float xx);
         void weightedMultiRegressionLinear(float** x,  float* y, float* weight, long nrItems,float* q,float* m, int nrPredictors);
         void weightedMultiRegressionLinear(const std::vector <std::vector <float>> &x, std::vector <float> &y, const std::vector <float> &weight, long nrItems,float* q,std::vector <float> &m, int nrPredictors);
+        void weightedMultiRegressionLinearWithoutConstantTerm(const std::vector <std::vector <float>> &x, std::vector <float> &y, const std::vector <float> &weight, long nrItems,std::vector <float> &m, int nrPredictors);
         void multiRegressionLinear(float** x,  float* y, long nrItems,float* q,float* m, int nrPredictors);
         void linearRegression(float* x, float* y, long nrItems, bool zeroIntercept, float* y_intercept, float* mySlope, float* r2);
         void linearRegression( std::vector<float> x,  std::vector<float> y, long nrItems, bool zeroIntercept, float* y_intercept, float* mySlope, float* r2);