linear regression Octave and C code

git-svn-id: https://svn.code.sf.net/p/freetel/code@2111 01035d8c-6547-0410-b346-abe4f91aad63

diff --git a/codec2-dev/octave/tlinreg.m b/codec2-dev/octave/tlinreg.m
new file mode 100644 (file)
index 0000000..81656e6
--- /dev/null
+++ b/codec2-dev/octave/tlinreg.m
@@ -0,0 +1,40 @@
+
+1;
+
+function [m b] = linreg(x,y,n)
+  sumx = 0.0;   % sum of x
+  sumx2 = 0.0;  % sum of x^2
+  sumxy = 0.0;  % sum of x * y
+  sumy = 0.0;   % sum of y
+  sumy2 = 0.0;  % sum of y**2
+
+  for i=1:n   
+    sumx  += x(i);       
+    sumx2 += x(i)^2;  
+    sumxy += x(i) * y(i);
+    sumy  += y(i);      
+    sumy2 += y(i)^2; 
+  end 
+
+  denom = (n * sumx2 - sumx*sumx);
+
+  if denom == 0
+    % singular matrix. can't solve the problem.
+    m = 0;
+    b = 0;
+  else
+    m = (n * sumxy  -  sumx * sumy) / denom;
+    b = (sumy * sumx2  -  sumx * sumxy) / denom;
+  end
+
+endfunction
+
+x = [1 2 7 8];
+y = [ -0.70702 + 0.70708i 0.77314 - 0.63442i -0.98083 + 0.19511i 0.99508 - 0.09799i] .* [1 -1 1 -1];
+[m b] = linreg(x,y,4);
+ 
+x = [1 2 3 4 5 6 7 8];
+yFit = m * x + b;
+figure(1)
+plot(yFit,'r*')
+axis([-1.5 1.5 -1.5 1.5])

diff --git a/codec2-dev/src/linreg.c b/codec2-dev/src/linreg.c
new file mode 100644 (file)
index 0000000..e91a69c
--- /dev/null
+++ b/codec2-dev/src/linreg.c
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: linreg.c
+  AUTHOR......: David Rowe
+  DATE CREATED: April 2015
+
+  Linear regression C module based on:
+
+    http://stackoverflow.com/questions/5083465/fast-efficient-least-squares-fit-algorithm-in-c
+
+  Use:
+
+    $ gcc linreg.c -o linreg -D__UNITTEST__ -Wall
+    $ ./linreg
+
+    Then compare yfit results with octave/tlinreg.m
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2015 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "linreg.h"
+#include "comp_prim.h"
+
+
+void linreg(COMP *m, COMP *b, float x[], COMP y[], int n)
+{
+    float  sumx  = 0.0;         /* sum of x          */
+    float  sumx2 = 0.0;         /* sum of x^2        */
+    COMP   sumxy = {0.0,0.0};   /* sum of x * y      */
+    COMP   sumy  = {0.0,0.0};   /* sum of y          */
+    COMP   sumy2 = {0.0,0.0};   /* sum of y**2       */
+    float  denom;
+    COMP   zero;
+    int    i;
+
+    for (i=0; i<n; i++) {   
+        sumx  += x[i];       
+        sumx2 += x[i]*x[i];  
+        sumxy = cadd(sumxy, fcmult(x[i], y[i]));
+        sumy  = cadd(sumy, y[i]);      
+        sumy2 = cadd(sumy2, cmult(y[i],y[i])); 
+    }
+
+  denom = (n * sumx2 - sumx*sumx);
+
+  if (denom == 0) {
+      /* singular matrix. can't solve the problem */
+      zero.real = 0.0; zero.imag = 0.0;
+      *m = zero;
+      *b = zero;
+  } else {
+      *m = fcmult(1.0/denom, cadd(fcmult(n, sumxy), cneg(fcmult(sumx,sumy))));
+      *b = fcmult(1.0/denom, cadd(fcmult(sumx2,sumy), cneg(fcmult(sumx, sumxy))));
+  }
+}
+
+
+#ifdef __UNITTEST__
+
+
+static float x[] = {1, 2, 7, 8};
+
+static COMP  y[] = {
+  {-0.70702,  0.70708},
+  {-0.77314,  0.63442},
+  {-0.98083,  0.19511},
+  {-0.99508,  0.09799}
+};
+
+
+int main(void) {
+    float  x1;
+    COMP   m,b,yfit;
+  
+    linreg(&m, &b, x, y, sizeof(x)/sizeof(float));
+
+    for (x1=1; x1<=8; x1++) {
+        yfit = cadd(fcmult(x1, m),b);
+        printf("%f %f\n", yfit.real, yfit.imag); 
+    }
+
+    return 0;
+}
+
+#endif
+

diff --git a/codec2-dev/src/linreg.h b/codec2-dev/src/linreg.h
new file mode 100644 (file)
index 0000000..c86dfad
--- /dev/null
+++ b/codec2-dev/src/linreg.h
@@ -0,0 +1,35 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: linreg.h
+  AUTHOR......: David Rowe
+  DATE CREATED: April 2015
+
+  Linear regression C module based
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2015 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __LINREG__
+#define __LINREG__
+
+#include "comp.h"
+
+void linreg(COMP *m, COMP *b, float x[], COMP y[], int n);
+
+#endif