Полиномиальная сила подгонки от первой степени до нуля

я нашел хороший код для выполнения полиномиальной подгонки методом наименьших квадратов на основе GSL.

я использую его с 3 степенями: y = Cx² + Bx + A.

В моем приложении я знаю, что A должен быть равен нулю. Можно ли изменить алгоритм так, чтобы A всегда было равно нулю?

bool polynomialfit(int obs, int degree, 
           double *dx, double *dy, double *store) /* n, p */
{
  gsl_multifit_linear_workspace *ws;
  gsl_matrix *cov, *X;
  gsl_vector *y, *c;
  double chisq;

  int i, j;

  X = gsl_matrix_alloc(obs, degree);
  y = gsl_vector_alloc(obs);
  c = gsl_vector_alloc(degree);
  cov = gsl_matrix_alloc(degree, degree);

  for(i=0; i < obs; i++) {
    gsl_matrix_set(X, i, 0, 1.0);
    for(j=0; j < degree; j++) {
      gsl_matrix_set(X, i, j, pow(dx[i], j));
    }
    gsl_vector_set(y, i, dy[i]);
  }

  ws = gsl_multifit_linear_alloc(obs, degree);
  gsl_multifit_linear(X, y, c, cov, &chisq, ws);

  /* store result ... */
  for(i=0; i < degree; i++)
  {
    store[i] = gsl_vector_get(c, i);
  }

  gsl_multifit_linear_free(ws);
  gsl_matrix_free(X);
  gsl_matrix_free(cov);
  gsl_vector_free(y);
  gsl_vector_free(c);
  return true; /* we do not "analyse" the result (cov matrix mainly)
  to know if the fit is "good" */
}

c++ math least-squares curve-fitting gsl
person Benjamin    schedule 16.11.2013    source источник
comment
Если ваши данные совместимы с предыдущим предположением, то минимизация действительно подтвердит это. Нет необходимости устанавливать A = 0. Но если вы хотите сравнить совпадения, GSL ясно утверждает, что X представляет собой матрицу n на p, где p = количество неизвестных параметров. Затем вам просто нужно удалить строку, связанную с x ^ 0, чтобы установить A = 0   -  person Vivian Miranda    schedule 16.11.2013

Ответы (2)


arrow_upward
0
arrow_downward

Вы можете заменить y на y' = y/x, а затем выполнить подбор полинома 1-й степени y'= Cx + B?

(если в вашем наборе данных присутствует точка x = 0, вы должны удалить ее, но эта точка не улучшает соответствие, если вы хотите применить ограничение A = 0, вы все равно можете использовать его для повторного вычисления согласия)

person Paweł Kordowski    schedule 16.11.2013

arrow_upward
0
arrow_downward

В коде, который вы разместили, есть этот цикл:

for(j=0; j < degree; j++) {
   gsl_matrix_set(X, i, j, pow(dx[i], j));
}

и функция pow вычисляет термины x^j, вы должны «игнорировать» термин, где j==0.

У меня нет доступа к GSL, поэтому следующее просто пришло мне в голову и не проверено:

bool polynomialfit(int obs, int polynom_degree, 
           double *dx, double *dy, double *store) /* n, p */
{
  gsl_multifit_linear_workspace *ws;
  gsl_matrix *cov, *X;
  gsl_vector *y, *c;
  double chisq;

  int i, j;
  int degree = polynom_degree - 1;

  X = gsl_matrix_alloc(obs, degree);
  y = gsl_vector_alloc(obs);
  c = gsl_vector_alloc(degree);
  cov = gsl_matrix_alloc(degree, degree);

  for(i=0; i < obs; i++) {
    gsl_matrix_set(X, i, 0, 1.0);
    for(j=0; j < degree; j++) {
      gsl_matrix_set(X, i, j, pow(dx[i], j+1));
    }
    gsl_vector_set(y, i, dy[i]);
  }

  ws = gsl_multifit_linear_alloc(obs, degree);
  gsl_multifit_linear(X, y, c, cov, &chisq, ws);

  /* store result ... */
  for(i=0; i < degree; i++)
  {
    store[i] = gsl_vector_get(c, i);
  }

  gsl_multifit_linear_free(ws);
  gsl_matrix_free(X);
  gsl_matrix_free(cov);
  gsl_vector_free(y);
  gsl_vector_free(c);
  return true; /* we do not "analyse" the result (cov matrix mainly)
  to know if the fit is "good" */
}

Чтобы соответствовать y=c*x*x+b*x, вы должны вызвать его с polynom_degree, установленным на 3.

Вы также можете ознакомиться с теорией:

Вайсштейн, Эрик В. «Подгонка методом наименьших квадратов - многочлен». Из MathWorld — веб-ресурса Wolfram. http://mathworld.wolfram.com/LeastSquaresFittingPolynomial.html

person Alessandro Jacopson    schedule 16.11.2013