Add argument specifying column containing dependent variable.

diff --git a/lib/linreg/sweep.c b/lib/linreg/sweep.c
index 000a5e3e0f740ed6b26957ccb2ee5703649db1e6..7d7dfc917fb3d846bdb29d7cc34564a8ea06b1d5 100644 (file)
--- a/lib/linreg/sweep.c
+++ b/lib/linreg/sweep.c
@@ -60,58 +60,78 @@
 
    The matrix A is assumed to be symmetric, so the sweep operation is
    performed only for the upper triangle of A.
+
+   LAST_COL is considered to be the final column in the augmented matrix,
+   that is, the column to the right of the '=' sign of the system.
  */
 
 int
-reg_sweep (gsl_matrix * A)
+reg_sweep (gsl_matrix * A, int last_col)
 {
   double sweep_element;
   double tmp;
   int i;
   int j;
   int k;
+  int row_i;
+  int row_k;
+  int row_j;
+  int *ordered_cols;
   gsl_matrix *B;
 
   if (A != NULL)
     {
       if (A->size1 == A->size2)
        {
+         ordered_cols = malloc (A->size1 * sizeof (*ordered_cols));
+         for (i = 0; i < last_col; i++)
+           {
+             ordered_cols[i] = i;
+           }
+         for (i = last_col + 1; i < A->size1 - 1; i++)
+           {
+             ordered_cols[i - 1] = i;
+           }
+         ordered_cols[A->size1 - 1] = last_col;
          B = gsl_matrix_alloc (A->size1, A->size2);
          for (k = 0; k < (A->size1 - 1); k++)
            {
-             sweep_element = gsl_matrix_get (A, k, k);
+             row_k = ordered_cols[k];
+             sweep_element = gsl_matrix_get (A, row_k, row_k);
              if (fabs (sweep_element) > GSL_DBL_MIN)
                {
                  tmp = -1.0 / sweep_element;
-                 gsl_matrix_set (B, k, k, tmp);
+                 gsl_matrix_set (B, row_k, row_k, tmp);
                  /*
                     Rows before current row k.
                   */
                  for (i = 0; i < k; i++)
                    {
+                     row_i = ordered_cols[i];
                      for (j = i; j < A->size2; j++)
                        {
+                         row_j = ordered_cols[j];
                          /*
                             Use only the upper triangle of A.
                           */
-                         if (j < k)
+                         if (row_j < row_k)
                            {
-                             tmp = gsl_matrix_get (A, i, j) -
-                               gsl_matrix_get (A, i, k)
-                               * gsl_matrix_get (A, j, k) / sweep_element;
-                             gsl_matrix_set (B, i, j, tmp);
+                             tmp = gsl_matrix_get (A, row_i, row_j) -
+                               gsl_matrix_get (A, row_i, row_k)
+                               * gsl_matrix_get (A, row_j, row_k) / sweep_element;
+                             gsl_matrix_set (B, row_i, row_j, tmp);
                            }
-                         else if (j > k)
+                         else if (row_j > row_k)
                            {
-                             tmp = gsl_matrix_get (A, i, j) -
-                               gsl_matrix_get (A, i, k)
-                               * gsl_matrix_get (A, k, j) / sweep_element;
-                             gsl_matrix_set (B, i, j, tmp);
+                             tmp = gsl_matrix_get (A, row_i, row_j) -
+                               gsl_matrix_get (A, row_i, row_k)
+                               * gsl_matrix_get (A, row_k, row_j) / sweep_element;
+                             gsl_matrix_set (B, row_i, row_j, tmp);
                            }
                          else
                            {
-                             tmp = gsl_matrix_get (A, i, k) / sweep_element;
-                             gsl_matrix_set (B, i, j, tmp);
+                             tmp = gsl_matrix_get (A, row_i, row_k) / sweep_element;
+                             gsl_matrix_set (B, row_i, row_j, tmp);
                            }
                        }
                    }
@@ -120,31 +140,37 @@ reg_sweep (gsl_matrix * A)
                   */
                  for (j = k + 1; j < A->size1; j++)
                    {
-                     tmp = gsl_matrix_get (A, k, j) / sweep_element;
-                     gsl_matrix_set (B, k, j, tmp);
+                     row_j = ordered_cols[j];
+                     tmp = gsl_matrix_get (A, row_k, row_j) / sweep_element;
+                     gsl_matrix_set (B, row_k, row_j, tmp);
                    }
                  /*
                     Rows after the current row k.
                   */
                  for (i = k + 1; i < A->size1; i++)
                    {
+                     row_i = ordered_cols[i];
                      for (j = i; j < A->size2; j++)
                        {
-                         tmp = gsl_matrix_get (A, i, j) -
-                           gsl_matrix_get (A, k, i)
-                           * gsl_matrix_get (A, k, j) / sweep_element;
-                         gsl_matrix_set (B, i, j, tmp);
+                         row_j = ordered_cols[j];
+                         tmp = gsl_matrix_get (A, row_i, row_j) -
+                           gsl_matrix_get (A, row_k, row_i)
+                           * gsl_matrix_get (A, row_k, row_j) / sweep_element;
+                         gsl_matrix_set (B, row_i, row_j, tmp);
                        }
                    }
                }
              for (i = 0; i < A->size1; i++)
                for (j = i; j < A->size2; j++)
                  {
-                   gsl_matrix_set (A, i, j, gsl_matrix_get (B, i, j));
+                   row_i = ordered_cols[i];
+                   row_j = ordered_cols[j];
+                   gsl_matrix_set (A, row_i, row_j, gsl_matrix_get (B, row_i, row_j));
                  }
            }
          gsl_matrix_free (B);
          return GSL_SUCCESS;
+         free (ordered_cols);
        }
       return GSL_ENOTSQR;
     }

diff --git a/lib/linreg/sweep.h b/lib/linreg/sweep.h
index 7bdd5c7c6c2a5f1f7dc43be81ecdb0c74db95f24..a1567a166222104f6ad839b2670c1d120f355f73 100644 (file)
--- a/lib/linreg/sweep.h
+++ b/lib/linreg/sweep.h
@@ -65,6 +65,6 @@
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_math.h>
 
-int reg_sweep (gsl_matrix * A);
+int reg_sweep (gsl_matrix *, int);
 
 #endif

diff --git a/src/math/linreg.c b/src/math/linreg.c
index e74646b97145b44720a2419b3bf5ff0a1f64cc1d..fe01ac39e4941969b8664d0dd8d08f7cf7c57d27 100644 (file)
--- a/src/math/linreg.c
+++ b/src/math/linreg.c
@@ -74,6 +74,7 @@ linreg_alloc (const struct variable *depvar, const struct variable **indep_vars,
   c = xmalloc (sizeof (*c));
   c->depvar = depvar;
   c->indep_vars = xnmalloc (p, sizeof (*indep_vars));
+  c->dependent_column = p;
   for (i = 0; i < p; i++)
     {
       c->indep_vars[i] = indep_vars[i];
@@ -375,7 +376,7 @@ linreg_fit (const gsl_matrix *cov, linreg *l)
       gsl_matrix *params;
       params = gsl_matrix_calloc (cov->size1, cov->size2);
       gsl_matrix_memcpy (params, cov);
-      reg_sweep (params);
+      reg_sweep (params, l->dependent_column);
       post_sweep_computations (l, params);  
       gsl_matrix_free (params);
     }
@@ -460,3 +461,9 @@ linreg_get_depvar_mean (linreg *c)
 {
   return c->depvar_mean;
 }
+static
+linreg_set_dependent_column (linreg *c, int column)
+{
+  c->dependent_column = column;
+}
+

diff --git a/src/math/linreg.h b/src/math/linreg.h
index 349d5a909c300c4713ee14a0196d967abbbcb5b6..31bf5651051564774d0417571a99f929a562f270 100644 (file)
--- a/src/math/linreg.h
+++ b/src/math/linreg.h
@@ -142,6 +142,7 @@ struct linreg_struct
 
   struct variable *pred;
   struct variable *resid;
+  int dependent_column; /* Column containing the dependent variable. Defaults to last column. */
 };
 
 typedef struct linreg_struct linreg;