/* * regression.c: Statistical regression functions. * * Authors: * Morten Welinder * Andrew Chatham * Daniel Carrera */ #include #include "gnumeric.h" #include "regression.h" #include "rangefunc.h" #include "mathfunc.h" #include #include #include #include #undef DEBUG_NEAR_SINGULAR #define ALLOC_MATRIX(var,dim1,dim2) \ do { int _i, _d1, _d2; \ _d1 = (dim1); \ _d2 = (dim2); \ (var) = g_new (gnm_float *, _d1); \ for (_i = 0; _i < _d1; _i++) \ (var)[_i] = g_new (gnm_float, _d2); \ } while (0) #define FREE_MATRIX(var,dim1,dim2) \ do { int _i, _d1; \ _d1 = (dim1); \ for (_i = 0; _i < _d1; _i++) \ g_free ((var)[_i]); \ g_free (var); \ } while (0) #define COPY_MATRIX(dst,src,dim1,dim2) \ do { int _i, _j, _d1, _d2; \ _d1 = (dim1); \ _d2 = (dim2); \ for (_i = 0; _i < _d1; _i++) \ for (_j = 0; _j < _d2; _j++) \ (dst)[_i][_j] = (src)[_i][_j]; \ } while (0) #define PRINT_MATRIX(var,dim1,dim2) \ do { \ int _i, _j, _d1, _d2; \ _d1 = (dim1); \ _d2 = (dim2); \ for (_i = 0; _i < _d1; _i++) \ { \ for (_j = 0; _j < _d2; _j++) \ fprintf (stderr, " %19.10" GNUM_FORMAT_g, (var)[_i][_j]); \ fprintf (stderr, "\n"); \ } \ } while (0) /* * ---> j * * | ******** * | ******** * | ******** A[i][j] * v ******** * ******** * i ******** * ******** * ******** * */ /* ------------------------------------------------------------------------- */ /* Returns in res the solution to the equation L * U * res = P * b. This function is adapted from pseudocode in Introduction to Algorithms_. Cormen, Leiserson, and Rivest. p. 753. MIT Press, 1990. */ static void backsolve (gnm_float **LU, int *P, gnm_float *b, int n, gnm_float *res) { int i, j; for (i = 0; i < n; i++) { res[i] = b[P[i]]; for (j = 0; j < i; j++) res[i] -= LU[i][j] * res[j]; } for (i = n - 1; i >= 0; i--) { for (j = i + 1; j < n; j++) res[i] -= LU[i][j] * res[j]; res[i] /= LU[i][i]; } } static RegressionResult rescale (gnm_float **A, gnm_float *b, int n, gnm_float *pdet) { int i; *pdet = 1; for (i = 0; i < n; i++) { int j, expn; gnm_float scale, max; (void)range_maxabs (A[i], n, &max); if (max == 0) return REG_singular; /* Use a power of 2 near sqrt (max) as scale. */ (void)frexpgnum (sqrtgnum (max), &expn); scale = ldexpgnum (1, expn); #ifdef DEBUG_NEAR_SINGULAR printf ("scale[%d]=%" GNUM_FORMAT_g "\n", i, scale); #endif *pdet *= scale; b[i] /= scale; for (j = 0; j < n; j++) A[i][j] /= scale; } return REG_ok; } /* * Performs an LUP Decomposition; LU and P must already be allocated. * A is not destroyed. * * This function is adapted from pseudocode in * _Introduction to Algorithms_. Cormen, Leiserson, and Rivest. * p 759. MIT Press, 1990. * * A rescaling of rows is done and the b_scaled vector is scaled * accordingly. */ static RegressionResult LUPDecomp (gnm_float **A, gnm_float **LU, int *P, int n, gnm_float *b_scaled, gnm_float *pdet) { int i, j, k, tempint; gnm_float highest = 0; gnm_float lowest = GNUM_MAX; gnm_float cond; gboolean odd_parity = FALSE; gnm_float det = 1; COPY_MATRIX (LU, A, n, n); for (j = 0; j < n; j++) P[j] = j; *pdet = 0; #ifdef DEBUG_NEAR_SINGULAR PRINT_MATRIX (LU, n, n); #endif { RegressionResult err = rescale (LU, b_scaled, n, &det); if (err != REG_ok) return err; } for (i = 0; i < n; i++) { gnm_float max = 0; int mov = -1; for (j = i; j < n; j++) if (gnumabs (LU[j][i]) > max) { max = gnumabs (LU[j][i]); mov = j; } #ifdef DEBUG_NEAR_SINGULAR PRINT_MATRIX (LU, n, n); printf ("max[%d]=%" GNUM_FORMAT_g " at %d\n", i, max, mov); #endif if (max == 0) return REG_singular; if (max > highest) highest = max; if (max < lowest) lowest = max; if (i != mov) { /*swap the two rows */ odd_parity = !odd_parity; tempint = P[i]; P[i] = P[mov]; P[mov] = tempint; for (j = 0; j < n; j++) { gnm_float temp = LU[i][j]; LU[i][j] = LU[mov][j]; LU[mov][j] = temp; } } for (j = i + 1; j < n; j++) { LU[j][i] /= LU[i][i]; for (k = i + 1; k < n; k++) LU[j][k] -= LU[j][i] * LU[i][k]; } } /* Calculate the determinant. */ if (odd_parity) det = -det; for (i = 0; i < n; i++) det *= LU[i][i]; *pdet = det; cond = (loggnum (highest) - loggnum (lowest)) / loggnum (2); #ifdef DEBUG_NEAR_SINGULAR printf ("cond=%.20" GNUM_FORMAT_g "\n", cond); #endif /* FIXME: make some science out of this. */ if (cond > GNUM_MANT_DIG * 0.75) return REG_near_singular_bad; else if (cond > GNUM_MANT_DIG * 0.50) return REG_near_singular_good; else return REG_ok; } static RegressionResult linear_solve (gnm_float **A, gnm_float *b, int n, gnm_float *res) { RegressionResult err; gnm_float **LU, *b_scaled; int *P; gnm_float det; if (n < 1) return REG_not_enough_data; /* Special case. */ if (n == 1) { gnm_float d = A[0][0]; if (d == 0) return REG_singular; res[0] = b[0] / d; return REG_ok; } /* Special case. */ if (n == 2) { gnm_float d = matrix_determinant (A, n); if (d == 0) return REG_singular; res[0] = (A[1][1] * b[0] - A[1][0] * b[1]) / d; res[1] = (A[0][0] * b[1] - A[0][1] * b[0]) / d; return REG_ok; } /* * Otherwise, use LUP-decomposition to find res such that * A res = b */ ALLOC_MATRIX (LU, n, n); P = g_new (int, n); b_scaled = g_new (gnm_float, n); memcpy (b_scaled, b, n * sizeof (gnm_float)); err = LUPDecomp (A, LU, P, n, b_scaled, &det); if (err == REG_ok || err == REG_near_singular_good) backsolve (LU, P, b_scaled, n, res); FREE_MATRIX (LU, n, n); g_free (P); g_free (b_scaled); return err; } gboolean matrix_invert (gnm_float **A, int n) { RegressionResult err; gnm_float **LU, *b_scaled, det; int *P; int i; gboolean res; if (n < 1) return FALSE; /* * Otherwise, use LUP-decomposition to find res such that * A res = b */ ALLOC_MATRIX (LU, n, n); P = g_new (int, n); b_scaled = g_new (gnm_float, n); for (i = 0; i < n; i++) b_scaled[i] = 1; err = LUPDecomp (A, LU, P, n, b_scaled, &det); if (err == REG_ok || err == REG_near_singular_good) { int i, j; gnm_float *b = g_new (gnm_float, n); gnm_float *w = g_new (gnm_float, n); for (i = 0; i < n; i++) { memset (b, 0, sizeof (gnm_float) * n); b[i] = b_scaled[i]; backsolve (LU, P, b, n, w); for (j = 0; j < n; j++) A[j][i] = w[j]; } g_free (w); g_free (b); res = TRUE; } else res = FALSE; FREE_MATRIX (LU, n, n); g_free (P); g_free (b_scaled); return res; } gnm_float matrix_determinant (gnm_float **A, int n) { RegressionResult err; gnm_float **LU, *b_scaled, det; int *P; if (n < 1) return 0; /* Special case. */ if (n == 1) return A[0][0]; /* Special case. */ if (n == 2) return A[0][0] * A[1][1] - A[1][0] * A[0][1]; /* * Otherwise, use LUP-decomposition to find res such that * A res = b */ ALLOC_MATRIX (LU, n, n); P = g_new (int, n); b_scaled = g_new0 (gnm_float, n); err = LUPDecomp (A, LU, P, n, b_scaled, &det); FREE_MATRIX (LU, n, n); g_free (P); g_free (b_scaled); return det; } /* ------------------------------------------------------------------------- */ static RegressionResult general_linear_regression (gnm_float **xss, int xdim, const gnm_float *ys, int n, gnm_float *result, regression_stat_t *regression_stat, gboolean affine) { gnm_float *xTy, **xTx; int i,j; RegressionResult regerr; if (regression_stat) memset (regression_stat, 0, sizeof (regression_stat_t)); if (xdim > n) return REG_not_enough_data; xTy = g_new (gnm_float, xdim); for (i = 0; i < xdim; i++) { const gnm_float *xs = xss[i]; register gnm_float res = 0; int j; if (xs == NULL) /* NULL represents a 1-vector. */ for (j = 0; j < n; j++) res += ys[j]; else for (j = 0; j < n; j++) res += xs[j] * ys[j]; xTy[i] = res; } ALLOC_MATRIX (xTx, xdim, xdim); for (i = 0; i < xdim; i++) { const gnm_float *xs1 = xss[i]; int j; for (j = 0; j <= i; j++) { const gnm_float *xs2 = xss[j]; gnm_float res = 0; int k; if (xs1 == NULL && xs2 == NULL) res = n; else if (xs1 == NULL) for (k = 0; k < n; k++) res += xs2[k]; else if (xs2 == NULL) for (k = 0; k < n; k++) res += xs1[k]; else for (k = 0; k < n; k++) res += xs1[k] * xs2[k]; xTx[i][j] = xTx[j][i] = res; } } regerr = linear_solve (xTx, xTy, xdim, result); if (regression_stat && (regerr == REG_ok || regerr == REG_near_singular_good)) { RegressionResult err2; gnm_float *residuals = g_new (gnm_float, n); gnm_float **LU, *one_scaled, det; int *P; int err; /* This should not fail since n >= 1. */ err = range_average (ys, n, ®ression_stat->ybar); g_assert (err == 0); /* FIXME: we ought to have a devsq variant that does not recompute the mean. */ if (affine) err = range_devsq (ys, n, ®ression_stat->ss_total); else err = range_sumsq (ys, n, ®ression_stat->ss_total); g_assert (err == 0); regression_stat->xbar = g_new (gnm_float, n); for (i = 0; i < xdim; i++) { if (xss[i]) { int err = range_average (xss[i], n, ®ression_stat->xbar[i]); g_assert (err == 0); } else { regression_stat->xbar[i] = 1; } } for (i = 0; i < n; i++) { residuals[i] = 0; for (j = 0; j < xdim; j++) { if (xss[j]) residuals[i] += xss[j][i] * result[j]; else residuals[i] += result[j]; /* If NULL, constant factor */ } residuals[i] = ys[i] - residuals[i]; } err = range_sumsq (residuals, n, ®ression_stat->ss_resid); g_assert (err == 0); regression_stat->sqr_r = (regression_stat->ss_total == 0) ? 1 : 1 - regression_stat->ss_resid / regression_stat->ss_total; /* FIXME: we want to guard against division by zero. */ regression_stat->adj_sqr_r = 1 - regression_stat->ss_resid * (n - 1) / ((n - xdim) * regression_stat->ss_total); regression_stat->var = (n == xdim) ? 0 : regression_stat->ss_resid / (n - xdim); ALLOC_MATRIX (LU, xdim, xdim); one_scaled = g_new (gnm_float, xdim); for (i = 0; i < xdim; i++) one_scaled[i] = 1; P = g_new (int, xdim); err2 = LUPDecomp (xTx, LU, P, xdim, one_scaled, &det); regression_stat->se = g_new (gnm_float, xdim); if (err2 == REG_ok || err2 == REG_near_singular_good) { gnm_float *e = g_new (gnm_float, xdim); /* Elementary vector */ gnm_float *inv = g_new (gnm_float, xdim); for (i = 0; i < xdim; i++) e[i] = 0; for (i = 0; i < xdim; i++) { e[i] = one_scaled[i]; backsolve (LU, P, e, xdim, inv); if (inv[i] < 0) { /* * If this happens, something is really * wrong, numerically. */ regerr = REG_near_singular_bad; } regression_stat->se[i] = sqrtgnum (regression_stat->var * inv[i]); e[i] = 0; } g_free (e); g_free (inv); } else { /* * This can happen for xdim == 2 as linear_solve does * not use LUPDecomp in that case. */ regerr = err2; for (i = 0; i < xdim; i++) regression_stat->se[i] = 0; } FREE_MATRIX (LU, xdim, xdim); g_free (P); g_free (one_scaled); regression_stat->t = g_new (gnm_float, xdim); for (i = 0; i < xdim; i++) regression_stat->t[i] = (regression_stat->se[i] == 0) ? gnm_pinf : result[i] / regression_stat->se[i]; regression_stat->df_resid = n - xdim; regression_stat->df_reg = xdim - (affine ? 1 : 0); regression_stat->df_total = regression_stat->df_resid + regression_stat->df_reg; regression_stat->F = (regression_stat->sqr_r == 1) ? gnm_pinf : ((regression_stat->sqr_r / regression_stat->df_reg) / (1 - regression_stat->sqr_r) * regression_stat->df_resid); regression_stat->ss_reg = regression_stat->ss_total - regression_stat->ss_resid; regression_stat->se_y = sqrtgnum (regression_stat->ss_total / n); regression_stat->ms_reg = (regression_stat->df_reg == 0) ? 0 : regression_stat->ss_reg / regression_stat->df_reg; regression_stat->ms_resid = (regression_stat->df_resid == 0) ? 0 : regression_stat->ss_resid / regression_stat->df_resid; g_free (residuals); } FREE_MATRIX (xTx, xdim, xdim); g_free (xTy); return regerr; } /* ------------------------------------------------------------------------- */ typedef struct { gnm_float min_x; gnm_float max_x; gnm_float min_y; gnm_float max_y; gnm_float mean_y; } point_cloud_measure_type; /* Takes the current 'sign' (res[0]) and 'c' (res[3]) from the calling * function, transforms xs to ln(sign*(x-c)), performs a simple * linear regression to find the best fitting 'a' (res[1]) and 'b' * (res[2]) for ys and transformed xs, and computes the sum of squared * residuals. * Needs 'sign' (i.e. +1 or -1) and 'c' so adjusted that (sign*(x-c)) is * positive for all xs. n must be > 0. These conditions are trusted to be * checked by the calling functions. * Is called often, so do not make it too slow. */ static int transform_x_and_linear_regression_log_fitting (gnm_float *xs, gnm_float *transf_xs, const gnm_float *ys, int n, gnm_float *res, point_cloud_measure_type *point_cloud) { int i; int result = REG_ok; gnm_float mean_transf_x, diff_x, resid_y; gnm_float sum1 = 0; gnm_float sum2 = 0; /* log (always > 0) */ for (i=0; imean_y); sum2 += diff_x * diff_x; } res[2] = sum1 / sum2; res[1] = point_cloud->mean_y - (res[2] * mean_transf_x); res[4] = 0; for (i=0; imax_x) - (point_cloud->min_x); /* Not needed here, but allocate it once for all subfunction calls */ transf_xs = g_new (gnm_float, n); /* Choose final accuracy of c with respect to range of xs. * Make accuracy be a whole power of 10. */ c_accuracy = log10gnum (x_range); if (c_accuracy < 0) if (modfgnum (c_accuracy, &c_accuracy_int) != 0) c_accuracy--; modfgnum (c_accuracy, &c_accuracy_int); c_accuracy = c_accuracy_int; c_accuracy = powgnum (10, c_accuracy); c_accuracy *= LOGFIT_C_ACCURACY; /* Determine sign. Take a c which is ``much to small'' since the part * of the curve cutting the point cloud is almost not bent. * If making c still smaller does not make things still worse, * assume that we have to change the direction of curve bending * by changing sign. */ c_step = x_range * LOGFIT_C_STEP_FACTOR; c_range = x_range * LOGFIT_C_RANGE_FACTOR; res[0] = 1; /* sign */ res[3] = point_cloud->min_x - c_range; temp_res[0] = 1; temp_res[3] = res[3] - c_step; transform_x_and_linear_regression_log_fitting (xs, transf_xs, ys, n, res, point_cloud); transform_x_and_linear_regression_log_fitting (xs, transf_xs, ys, n, temp_res, point_cloud); if (temp_res[4] <= res[4]) sign_plus_ok = 0; /* check again with new sign */ res[0] = -1; /* sign */ res[3] = point_cloud->max_x + c_range; temp_res[0] = -1; temp_res[3] = res[3] + c_step; transform_x_and_linear_regression_log_fitting (xs, transf_xs, ys, n, res, point_cloud); transform_x_and_linear_regression_log_fitting (xs, transf_xs, ys, n, temp_res, point_cloud); if (temp_res[4] <= res[4]) sign_minus_ok = 0; /* If not exactly one of plus or minus works, give up. * This happens in point clouds which are very weakly bent. */ if (sign_plus_ok && !sign_minus_ok) res[0] = 1; else if (sign_minus_ok && !sign_plus_ok) res[0] = -1; else { result = REG_invalid_data; goto out; } /* Start of fitted c-range. Rounded to final accuracy of c. */ c_offset = (res[0] == 1) ? point_cloud->min_x : point_cloud->max_x; c_offset = c_accuracy * ((res[0] == 1) ? floorgnum (c_offset / c_accuracy) : ceilgnum (c_offset /c_accuracy)); /* Now the adapting of c starts. Find a local minimum of sum * of squared residuals. */ /* First, catch some unsuitably shaped point clouds. */ res[3] = c_offset - res[0] * c_accuracy; temp_res[3] = c_offset - res[0] * 2 * c_accuracy; temp_res[0] = res[0]; transform_x_and_linear_regression_log_fitting (xs, transf_xs, ys, n, res, point_cloud); transform_x_and_linear_regression_log_fitting (xs, transf_xs, ys, n, temp_res, point_cloud); if (temp_res[4] >= res[4]) { result = REG_invalid_data; goto out; } /* After the above check, any minimum reached will be NOT at * the start of c-range (c_offset - sign * c_accuracy) */ c_start = c_offset; c_end = c_start - res[0] * c_range; c_dist = res[0] * (c_start - c_end) / 2; res[3] = c_end + res[0] * c_dist; do { c_dist /= 2; transform_x_and_linear_regression_log_fitting (xs, transf_xs, ys, n, res, point_cloud); temp_res[3] = res[3] + res[0] * c_dist; transform_x_and_linear_regression_log_fitting (xs, transf_xs, ys, n, temp_res, point_cloud); if (temp_res[4] <= res[4]) memcpy (res, temp_res, 5 * sizeof (gnm_float)); else { temp_res[3] = res[3] - res[0] * c_dist; transform_x_and_linear_regression_log_fitting (xs, transf_xs, ys, n, temp_res, point_cloud); if (temp_res[4] <= res[4]) memcpy (res, temp_res, 5*sizeof (gnm_float)); } } while (c_dist > c_accuracy); res[3] = c_accuracy * gnumeric_fake_round (res[3] / c_accuracy); transform_x_and_linear_regression_log_fitting (xs, transf_xs, ys, n, res, point_cloud); if ((res[0] * (res[3] - c_end)) < (1.1 * c_accuracy)) { /* Allowing for some inaccuracy, we are at the end of the * range, so this is probably no local minimum. * The start of the range has been checked above. */ result = REG_invalid_data; goto out; } out: g_free (transf_xs); g_free (temp_res); return result; } /* ------------------------------------------------------------------------- */ /* Please refer to description in regression.h. */ RegressionResult linear_regression (gnm_float **xss, int dim, const gnm_float *ys, int n, gboolean affine, gnm_float *res, regression_stat_t *regression_stat) { RegressionResult result; g_return_val_if_fail (dim >= 1, REG_invalid_dimensions); g_return_val_if_fail (n >= 1, REG_invalid_dimensions); if (affine) { gnm_float **xss2; xss2 = g_new (gnm_float *, dim + 1); xss2[0] = NULL; /* Substitute for 1-vector. */ memcpy (xss2 + 1, xss, dim * sizeof (gnm_float *)); result = general_linear_regression (xss2, dim + 1, ys, n, res, regression_stat, affine); g_free (xss2); } else { res[0] = 0; result = general_linear_regression (xss, dim, ys, n, res + 1, regression_stat, affine); } return result; } /* ------------------------------------------------------------------------- */ /* Please refer to description in regression.h. */ RegressionResult exponential_regression (gnm_float **xss, int dim, const gnm_float *ys, int n, gboolean affine, gnm_float *res, regression_stat_t *regression_stat) { gnm_float *log_ys; RegressionResult result; int i; g_return_val_if_fail (dim >= 1, REG_invalid_dimensions); g_return_val_if_fail (n >= 1, REG_invalid_dimensions); log_ys = g_new (gnm_float, n); for (i = 0; i < n; i++) if (ys[i] > 0) log_ys[i] = loggnum (ys[i]); else { result = REG_invalid_data; goto out; } if (affine) { gnm_float **xss2; xss2 = g_new (gnm_float *, dim + 1); xss2[0] = NULL; /* Substitute for 1-vector. */ memcpy (xss2 + 1, xss, dim * sizeof (gnm_float *)); result = general_linear_regression (xss2, dim + 1, log_ys, n, res, regression_stat, affine); g_free (xss2); } else { res[0] = 0; result = general_linear_regression (xss, dim, log_ys, n, res + 1, regression_stat, affine); } if (result == 0) for (i = 0; i < dim + 1; i++) res[i] = expgnum (res[i]); out: g_free (log_ys); return result; } /* ------------------------------------------------------------------------- */ /* Please refer to description in regression.h. */ RegressionResult logarithmic_regression (gnm_float **xss, int dim, const gnm_float *ys, int n, gboolean affine, gnm_float *res, regression_stat_t *regression_stat) { gnm_float **log_xss; RegressionResult result; int i, j; g_return_val_if_fail (dim >= 1, REG_invalid_dimensions); g_return_val_if_fail (n >= 1, REG_invalid_dimensions); ALLOC_MATRIX (log_xss, dim, n); for (i = 0; i < dim; i++) for (j = 0; j < n; j++) if (xss[i][j] > 0) log_xss[i][j] = loggnum (xss[i][j]); else { result = REG_invalid_data; goto out; } if (affine) { gnm_float **log_xss2; log_xss2 = g_new (gnm_float *, dim + 1); log_xss2[0] = NULL; /* Substitute for 1-vector. */ memcpy (log_xss2 + 1, log_xss, dim * sizeof (gnm_float *)); result = general_linear_regression (log_xss2, dim + 1, ys, n, res, regression_stat, affine); g_free (log_xss2); } else { res[0] = 0; result = general_linear_regression (log_xss, dim, ys, n, res + 1, regression_stat, affine); } out: FREE_MATRIX (log_xss, dim, n); return result; } /* ------------------------------------------------------------------------- */ /* Please refer to description in regression.h. */ RegressionResult logarithmic_fit (gnm_float *xs, const gnm_float *ys, int n, gnm_float *res) { point_cloud_measure_type point_cloud_measures; int i, result; gboolean more_2_y = 0, more_2_x = 0; /* Store useful measures for using them here and in subfunctions. * The checking of n is paranoid -- the calling function should * have cared for that. */ g_return_val_if_fail (n > 2, REG_invalid_dimensions); result = range_min (xs, n, &(point_cloud_measures.min_x)); result = range_max (xs, n, &(point_cloud_measures.max_x)); result = range_min (ys, n, &(point_cloud_measures.min_y)); result = range_max (ys, n, &(point_cloud_measures.max_y)); result = range_average (ys, n, &(point_cloud_measures.mean_y)); /* Checking of error conditions. */ /* less than 2 different ys or less than 2 different xs */ g_return_val_if_fail (((point_cloud_measures.min_y != point_cloud_measures.max_y) && (point_cloud_measures.min_x != point_cloud_measures.max_x)), REG_invalid_data); /* less than 3 different ys */ for (i=0; ise = NULL; regression_stat->t = NULL; regression_stat->xbar = NULL; return regression_stat; } /* ------------------------------------------------------------------------- */ void regression_stat_destroy (regression_stat_t *regression_stat) { g_return_if_fail (regression_stat != NULL); if (regression_stat->se) g_free(regression_stat->se); if (regression_stat->t) g_free(regression_stat->t); if (regression_stat->xbar) g_free(regression_stat->xbar); g_free (regression_stat); } /* ------------------------------------------------------------------------- */ #define DELTA 0.01 /* FIXME: I pulled this number out of my hat. * I need some testing to pick a sensible value. */ #define MAX_STEPS 200 /* * SYNOPSIS: * result = derivative( f, &df, x, par, i) * * Approximates the partial derivative of a given function, at (x;params) * with respect to the parameter indicated by ith parameter. The resulst * is stored in 'df'. * * See the header file for more information. */ static RegressionResult derivative (RegressionFunction f, gnm_float *df, gnm_float *x, /* Only one point, not the whole data set. */ gnm_float *par, int index) { gnm_float y1, y2; RegressionResult result; gnm_float par_save = par[index]; par[index] = par_save - DELTA; result = (*f) (x, par, &y1); if (result != REG_ok) { par[index] = par_save; return result; } par[index] = par_save + DELTA; result = (*f) (x, par, &y2); if (result != REG_ok) { par[index] = par_save; return result; } #ifdef DEBUG printf ("y1 = %lf\n", y1); printf ("y2 = %lf\n", y2); printf ("DELTA = %lf\n",DELTA); #endif *df = (y2 - y1) / (2 * DELTA); par[index] = par_save; return REG_ok; } /* * SYNOPSIS: * result = chi_squared (f, xvals, par, yvals, sigmas, x_dim, &chisq) * * / y - f(x ; par) \ 2 * 2 | i i | * Chi == Sum ( | ------------------ | ) * \ sigma / * i * * sigmas -> Measurement errors in the dataset (along the y-axis). * NULL means "no errors available", so they are all set to 1. * * x_dim -> Number of data points. * * This value is not very meaningful without the sigmas. However, it is * still useful for the fit. */ static RegressionResult chi_squared (RegressionFunction f, gnm_float ** xvals, /* The entire data set. */ gnm_float *par, gnm_float *yvals, /* Ditto. */ gnm_float *sigmas, /* Ditto. */ int x_dim, /* Number of data points. */ gnm_float *chisq) /* Chi Squared */ { int i; RegressionResult result; gnm_float tmp, y; *chisq = 0; for (i = 0; i < x_dim; i++) { result = f (xvals[i], par, &y); if (result != REG_ok) return result; tmp = (yvals[i] - y ) / (sigmas ? sigmas[i] : 1); *chisq += tmp * tmp; } return REG_ok; } /* * SYNOPSIS: * result = chi_derivative (f, &dchi, xvals, par, i, yvals, * sigmas, x_dim) * * This is a simple adaptation of the derivative() function specific to * the Chi Squared. */ static RegressionResult chi_derivative (RegressionFunction f, gnm_float *dchi, gnm_float **xvals, /* The entire data set. */ gnm_float *par, int index, gnm_float *yvals, /* Ditto. */ gnm_float *sigmas, /* Ditto. */ int x_dim) { gnm_float y1, y2; RegressionResult result; gnm_float par_save = par[index]; par[index] = par_save - DELTA; result = chi_squared (f, xvals, par, yvals, sigmas, x_dim, &y1); if (result != REG_ok) { par[index] = par_save; return result; } par[index] = par_save + DELTA; result = chi_squared (f, xvals, par, yvals, sigmas, x_dim, &y2); if (result != REG_ok) { par[index] = par_save; return result; } #ifdef DEBUG printf ("y1 = %lf\n", y1); printf ("y2 = %lf\n", y2); printf ("DELTA = %lf\n", DELTA); #endif *dchi = (y2 - y1) / (2 * DELTA); par[index] = par_save; return REG_ok; } /* * SYNOPSIS: * result = coefficient_matrix (A, f, xvals, par, yvals, sigmas, * x_dim, p_dim, r) * * RETURNS: * The coefficient matrix of the LM method. * * DETAIS: * The coefficient matrix matrix is defined by * * N 1 df df * A = Sum ( ------- -- -- ( i == j ? 1 + r : 1 ) a) * ij k=1 sigma^2 dp dp * k i j * * A -> p_dim X p_dim coefficient matrix. MUST ALREADY BE ALLOCATED. * * sigmas -> Measurement errors in the dataset (along the y-axis). * NULL means "no errors available", so they are all set to 1. * * x_dim -> Number of data points. * * p_dim -> Number of parameters. * * r -> Positive constant. It's value is altered during the LM procedure. */ static RegressionResult coefficient_matrix (gnm_float **A, /* Output matrix. */ RegressionFunction f, gnm_float **xvals, /* The entire data set. */ gnm_float *par, gnm_float *yvals, /* Ditto. */ gnm_float *sigmas, /* Ditto. */ int x_dim, /* Number of data points. */ int p_dim, /* Number of parameters. */ gnm_float r) { int i, j, k; RegressionResult result; gnm_float df_i, df_j; gnm_float sum, sigma; /* Notice that the matrix is symetric. */ for (i = 0; i < p_dim; i++) { for (j = 0; j <= i; j++) { sum = 0; for (k = 0; k < x_dim; k++) { result = derivative (f, &df_i, xvals[k], par, i); if (result != REG_ok) return result; result = derivative (f, &df_j, xvals[k], par, j); if (result != REG_ok) return result; sigma = (sigmas ? sigmas[k] : 1); sum += (df_i * df_j) / (sigma * sigma) * (i == j ? 1 + r : 1) ; } A[i][j] = A[j][i] = sum; } } return REG_ok; } /* * SYNOPSIS: * result = parameter_errors (f, xvals, par, yvals, sigmas, * x_dim, p_dim, errors) * * Returns the errors associated with the parameters. * If an error is infinite, it is set to -1. * * sigmas -> Measurement errors in the dataset (along the y-axis). * NULL means "no errors available", so they are all set to 1. * * x_dim -> Number of data points. * * p_dim -> Number of parameters. * * errors -> MUST ALREADY BE ALLOCATED. */ /* FIXME: I am not happy with the behaviour with infinite errors. */ static RegressionResult parameter_errors (RegressionFunction f, gnm_float **xvals, /* The entire data set. */ gnm_float *par, gnm_float *yvals, /* Ditto. */ gnm_float *sigmas, /* Ditto. */ int x_dim, /* Number of data points. */ int p_dim, /* Number of parameters. */ gnm_float *errors) { RegressionResult result; gnm_float **A; int i; ALLOC_MATRIX (A, p_dim, p_dim); result = coefficient_matrix (A, f, xvals, par, yvals, sigmas, x_dim, p_dim, 0); if (result == REG_ok) { for (i = 0; i < p_dim; i++) /* FIXME: these were "[i][j]" which makes no sense. */ errors[i] = (A[i][i] != 0 ? 1 / sqrtgnum (A[i][i]) : -1); } FREE_MATRIX (A, p_dim, p_dim); return result; } /* * SYNOPSIS: * result = non_linear_regression (f, xvals, par, yvals, sigmas, * x_dim, p_dim, &chi, errors) * * Returns the results of the non-linear regression from the given initial * values. * The resulting parameters are placed back into 'par'. * * PARAMETERS: * * sigmas -> Measurement errors in the dataset (along the y-axis). * NULL means "no errors available", so they are all set to 1. * * x_dim -> Number of data points. * * p_dim -> Number of parameters. * * errors -> MUST ALREADY BE ALLOCATED. These are the approximated standard * deviation for each parameter. * * chi -> Chi Squared of the final result. This value is not very * meaningful without the sigmas. */ RegressionResult non_linear_regression (RegressionFunction f, gnm_float **xvals, /* The entire data set. */ gnm_float *par, gnm_float *yvals, /* Ditto. */ gnm_float *sigmas, /* Ditto. */ int x_dim, /* Number of data points. */ int p_dim, /* Number of parameters. */ gnm_float *chi, gnm_float *errors) { gnm_float r = 0.001; /* Pick a conservative initial value. */ gnm_float *b, **A; gnm_float *dpar; gnm_float *tmp_par; gnm_float chi_pre, chi_pos, dchi; RegressionResult result; int i, count; result = chi_squared (f, xvals, par, yvals, sigmas, x_dim, &chi_pre); if (result != REG_ok) return result; ALLOC_MATRIX (A, p_dim, p_dim); dpar = g_new (gnm_float, p_dim); tmp_par = g_new (gnm_float, p_dim); b = g_new (gnm_float, p_dim); #ifdef DEBUG printf ("Chi Squared : %lf", chi_pre); #endif for (count = 0; count < MAX_STEPS; count++) { for (i = 0; i < p_dim; i++) { /* * d Chi * b == ----- * k d p * k */ result = chi_derivative (f, &dchi, xvals, par, i, yvals, sigmas, x_dim); if (result != REG_ok) goto out; b[i] = - dchi; } result = coefficient_matrix (A, f, xvals, par, yvals, sigmas, x_dim, p_dim, r); if (result != REG_ok) goto out; result = linear_solve (A, b, p_dim, dpar); if (result != REG_ok) goto out; for(i = 0; i < p_dim; i++) tmp_par[i] = par[i] + dpar[i]; result = chi_squared (f, xvals, par, yvals, sigmas, x_dim, &chi_pos); if (result != REG_ok) goto out; #ifdef DEBUG printf ("Chi Squared : %lf", chi_pre); printf ("Chi Squared : %lf", chi_pos); printf ("r : %lf", r); #endif if (chi_pos <= chi_pre + DELTA / 2) { /* There is improvement */ r /= 10; par = tmp_par; if (gnumabs (chi_pos - chi_pre) < DELTA) break; chi_pre = chi_pos; } else { r *= 10; } } result = parameter_errors (f, xvals, par, yvals, sigmas, x_dim, p_dim, errors); if (result != REG_ok) goto out; *chi = chi_pos; out: FREE_MATRIX (A, p_dim, p_dim); g_free (dpar); g_free (tmp_par); g_free (b); return result; }