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

я нашел хороший код для подбора полиномиальных наименьших квадратов на основе 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" */
}

1

Решение

Вы можете заменить y на y ‘= y / x, а затем выполнить подгонку полинома 1. степени y’ = Cx + B?

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

0

Другие решения

В размещенном вами коде есть этот цикл:

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

0

По вопросам рекламы [email protected]