(免责声明:我的数学很糟糕，而且我来自 JavaScript，所以对于任何不准确的地方，我深表歉意，并将尽力纠正它们。)
Rosetta Code 上的 example 展示了如何使用 gsl 计算系数。这是代码:
polifitgsl.h:
#ifndef _POLIFITGSL_H
#define _POLIFITGSL_H
#include <gsl/gsl_multifit.h>
#include <stdbool.h>
#include <math.h>
bool polynomialfit(int obs, int degree,
double *dx, double *dy, double *store); /* n, p */
#endif
polifitgsl.cpp:
#include "polifitgsl.h"
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++) {
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" */
}
main.cpp(注意我已经用我自己的替换了 x 和 y 的样本编号):
#include <stdio.h>
#include "polifitgsl.h"
#define NP 11
double x[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19};
double y[] = {98.02, 98.01, 98.01, 98.02, 97.98, 97.97, 97.96, 97.94, 97.96, 97.96, 97.97, 97.97, 97.94, 97.94, 97.94, 97.92, 97.96, 97.9, 97.85, 97.9};
#define DEGREE 3
double coeff[DEGREE];
int main()
{
int i;
polynomialfit(NP, DEGREE, x, y, coeff);
for(i=0; i < DEGREE; i++) {
printf("%lf\n", coeff[i]);
}
return 0;
}
这是输出:
98.030909
-0.016182
0.000909
所以这给了我系数。但我真正想要的是实际的拟合点。在 JavaScript 中，我使用 regression package 来计算点数:
var regression = require('regression');
var calculateRegression = function(values, degree) {
var data = [];
var regressionOutput;
var valuesCount = values.length;
var i = 0;
// Format the data in a way the regression library expects.
for (i = 0; i < valuesCount; i++) {
data[i] = [i, values[i]];
}
// Use the library to calculate the regression.
regressionOutput = regression('polynomial', data, degree);
return regressionOutput;
};
var y = [98.02, 98.01, 98.01, 98.02, 97.98, 97.97, 97.96, 97.94, 97.96, 97.96, 97.97, 97.97, 97.94, 97.94, 97.94, 97.92, 97.96, 97.9, 97.85, 97.9];
console.log(calculateRegression(y, 3));
产生:
{ equation:
[ 98.02987916431594,
-0.017378390369880512,
0.0015748071645344357,
-0.00005721503635571101 ],
points:
[ [ 0, 98.02987916431594 ],
[ 1, 98.01401836607424 ],
[ 2, 98.00096389194348 ],
[ 3, 97.9903724517055 ],
[ 4, 97.98190075514219 ],
[ 5, 97.97520551203543 ],
[ 6, 97.96994343216707 ],
[ 7, 97.96577122531896 ],
[ 8, 97.96234560127297 ],
[ 9, 97.959323269811 ],
[ 10, 97.95636094071487 ],
[ 11, 97.95311532376647 ],
[ 12, 97.94924312874768 ],
[ 13, 97.94440106544033 ],
[ 14, 97.93824584362629 ],
[ 15, 97.93043417308745 ],
[ 16, 97.92062276360569 ],
[ 17, 97.90846832496283 ],
[ 18, 97.89362756694074 ],
[ 19, 97.87575719932133 ] ],
string: 'y = 0x^3 + 0x^2 + -0.02x + 98.03' }
(注意 JavaScript 中存在 float 问题，因此数字并不完全准确。)
points 这是我想要使用 gsl 生成的内容。我有办法做到这一点吗？
最佳答案
Chad，如果我理解您的需求，您只需使用通过 polynomialfit 函数找到的系数，根据多项式方程计算值。计算系数后，您可以使用以下等式找到任何 x(对于 DEGREE = 3)的值 y:
y = x^2*(coeff2) + x*(coeff1) + coeff0
或 C 语法
y = x*x*coeff[2] + x*coeff[1] + coeff[0];
您可以按如下方式修改您的main.cpp(您应该将main.cpp 重命名为main.c 和polyfitgsl .cpp 到 polyfitgsl.c -- 因为它们都是 C 源文件，而不是 C++)
#include <stdio.h>
#include "polifitgsl.h"
#define NP 20
double x[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
10, 11, 12, 13, 14, 15, 16, 17, 18, 19};
double y[] = { 98.02, 98.01, 98.01, 98.02, 97.98,
97.97, 97.96, 97.94, 97.96, 97.96,
97.97, 97.97, 97.94, 97.94, 97.94,
97.92, 97.96, 97.9, 97.85, 97.9 };
#define DEGREE 3
double coeff[DEGREE];
int main (void)
{
int i;
polynomialfit (NP, DEGREE, x, y, coeff);
printf ("\n polynomial coefficients:\n\n");
for (i = 0; i < DEGREE; i++) {
printf (" coeff[%d] : %11.7lf\n", i, coeff[i]);
}
putchar ('\n');
printf (" computed values:\n\n x y\n");
for (i = 0; i < (int)(sizeof x/sizeof *x); i++)
printf (" %2d, %11.7lf\n", i,
i*i*coeff[2] + i*coeff[1] + coeff[0]);
putchar ('\n');
return 0;
}
它提供了以下输出:
示例使用/输出
$ ./bin/main
polynomial coefficients:
coeff[0] : 98.0132468
coeff[1] : -0.0053003
coeff[2] : -0.0000558
computed values:
x y
0, 98.0132468
1, 98.0078906
2, 98.0024229
3, 97.9968435
4, 97.9911524
5, 97.9853497
6, 97.9794354
7, 97.9734094
8, 97.9672718
9, 97.9610226
10, 97.9546617
11, 97.9481891
12, 97.9416049
13, 97.9349091
14, 97.9281016
15, 97.9211825
16, 97.9141517
17, 97.9070093
18, 97.8997553
19, 97.8923896
这似乎是你所要求的。查看它，比较您的值(value)观，如果您需要其他帮助，请告诉我。
进一步清理
只是为了进一步清理代码，因为不需要全局声明 x、y 或 coeff，并使用 enum 定义常量 DEGREE 和 NP，更简洁的版本是:
#include <stdio.h>
#include "polifitgsl.h"
enum { DEGREE = 3, NP = 20 };
int main (void)
{
int i;
double x[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
10, 11, 12, 13, 14, 15, 16, 17, 18, 19};
double y[] = { 98.02, 98.01, 98.01, 98.02, 97.98,
97.97, 97.96, 97.94, 97.96, 97.96,
97.97, 97.97, 97.94, 97.94, 97.94,
97.92, 97.96, 97.9, 97.85, 97.9 };
double coeff[DEGREE] = {0};
polynomialfit (NP, DEGREE, x, y, coeff);
printf ("\n polynomial coefficients:\n\n");
for (i = 0; i < DEGREE; i++)
printf (" coeff[%d] : %11.7lf\n", i, coeff[i]);
putchar ('\n');
printf (" computed values:\n\n x y\n");
for (i = 0; i < (int)(sizeof x/sizeof *x); i++)
printf (" %2d, %11.7lf\n", i,
i*i*coeff[2] + i*coeff[1] + coeff[0]);
putchar ('\n');
return 0;
}
关于c - 如何使用 gsl 计算多项式回归数据点？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/36522882/