4

短：

GNUPlot 比我的 GSL 代码更适合我的数据。为什么？

短的：

我现在有点困惑，所以我的问题可能措辞不是特别好......随着我的理解提高，我会编辑这个。

这个问题的原标题是：“g++ Compiling code with either -o1 -o2 or -o3 flags and floating point precision”

我相信我的代码正遭受数值不稳定的困扰。

GNUPlot 比我的 GSL 代码更适合我的数据，这令人惊讶，因为我相信 GNUPlot 也使用 GSL 库？

长：

我编写了一些使用 GNU 科学库 (GSL) 的 C/C++ 代码。我的代码使非线性函数适合非线性数据集。执行此操作的算法可能对浮点运算发生的顺序高度敏感，这是由于数值不准确性的性质导致数值舍入误差的累积。

问题：“这些可能是由使用优化标志之一运行的效果引起的-o1，-o2还是-o3？”

部分答案：我关闭了所有-oN标志并重新编译了我的代码，我的结果可能会发生少量变化，即：delta_x / x ~= 1.0e-3. 与 GNUPlot 相比，拟合度仍然很差。

我适合的功能：

我提供这些是为了向您展示正在发生的数字工作。我怀疑其中一些容易出现数值错误。

的典型值Yi将在 至 的范围0.0内1.0。t通常在 到 的范围0.0内200.0。（但在该范围的前半部分拟合度很差。）

// Block A - Weighted Residuals

double t = time; // whatever that may be
double s = sigma[ix]; // these are the errors, tried 1.0 and 0.001 with no affect on parameter values obtained

double Yi = (std::sqrt(rho) * r0) / std::sqrt((rho - r0*r0) * std::exp(-2.0 * t / tau) + r0*r0); // y value - model

gsl_vector_set(f, ix, (Yi - y[ix])/sigma[ix]); // weighted residual

// Block B - Jacobian

double MEM_A = std::exp(-2.0 * t / tau); // Tried to do some optimization here
double MEM_B = 1.0 - MEM_A; // Perhaps this is causing a problem?
double MEM_C = rho - r0 * r0;
double MEM_D = MEM_C * MEM_A + r0*r0;
double MEM_E = std::pow(MEM_D, 1.5);

double df_drho = (std::pow(r0, 3.0) * MEM_B) / (2.0 * std::sqrt(rho) * MEM_E);
double df_dr0 = (std::pow(rho, 1.5) * MEM_A) / MEM_E;
double df_dtau = -1.0 * (std::sqrt(rho) * r0 * MEM_C * MEM_A * t) / (tau * tau * MEM_E);

gsl_matrix_set(J, ix, 0, df_drho / s);
gsl_matrix_set(J, ix, 1, df_dr0 / s);
gsl_matrix_set(J, ix, 2, df_dtau / s);

这是一个图表，不是很好吗？

好吧，这里有一张图表，它比我用语言更好地解释了这个问题。您可以忽略绿线，它仅显示在运行拟合算法之前给出的初始参数，该算法会更改这些参数。

GNUPlot 拟合结果：

RHOFIT = 0.086173236829715 +- 2.61304934752193e-05
 R0FIT = 0.00395856812689133 +- 2.08898744280108e-05
TAUFIT = 11.7694359189233 +- 0.016094629240588

// Not sure how GNUPlot calculates errors - they do not appear to be the regular errors computed from the off diagonal elements of the LM matrix after fitting. (If you know a little about the Levenberg–Marquardt algorithm.)

C++ GSL 拟合结果：

rho    = 0.08551510 +/- ...
r0     = 0.00507645 +/- ... // Nonsense errors due to "not-real" errors on data points
tau    = 12.99114719 +/- ...

仔细检查后，您会发现粉红色和蓝色线条并没有以相当大的余量相互重叠。粉红线是许多人所说的“合身”。相比之下，蓝线并不是特别好。

我已经尝试使误差条（尽管它们对于所有点的大小都相同——它们不是“真正的”误差条，只是人造的）更小——这没有帮助，只会改变每个点的卡方值和相关错误拟合后的参数。

图 1

进一步的随机想法：

  • 我建的 GSL 错了吗？
  • Gnuplot 将数据集拆分为小块，以使加在一起的数字保持大致相同的数量级？（有点像 FFT 的工作原理。）

GSL 拟合输出：

iter: 0 x = 0.1 0.001 10 |f(x)| = 12487.8
status = success
iter: 1 x = 0.0854247 0.00323946 13.2064 |f(x)| = 10476.9
dx vector: -0.0145753, 0.00223946, 3.20642
status = success
iter: 2 x = 0.0854309 0.00576809 13.7443 |f(x)| = 3670.4
dx vector: 6.18836e-06, 0.00252863, 0.537829
chisq/dof = 6746.03
rho    = 0.08543089 +/- 0.00013518
r0     = 0.00576809 +/- 0.00013165
tau    = 13.74425294 +/- 0.09012196
4

1 回答 1

2

我遇到了这个页面，因为我遇到了完全相同的问题。我需要用 GSL 拟合一个函数，以前没有这样做过，所以我将结果与 gnuplot 的拟合例程进行比较。在我的例子中，我正在为星系功率谱的一部分拟合一个简单的幂律，而 GSL 给我的拟合 chi^2/DoF 约为 6。

为了解决这个问题，我发现我很粗心，我的数据点的x值与正在评估拟合函数的x值不匹配。解决此问题的最简单方法是从数据值创建样条曲线，然后在将评估拟合函数的相同x值处评估样条曲线。例如：

#include <gsl/gsl_spline.h>
    .
    .
    .
std::vector< double > xvals;
std::vector< double > yvals;

fin.open("SomeDataFile.dat", std::ios::in);
while (!fin.eof()) {
    double x, y;
    fin >> x >> y;
    xvals.push_back(x);
    yvals.push_back(y);
}
fin.close();

gsl_spline *Y = gsl_spline_alloc(gsl_interp_cspline, yvals.size());
gsl_interp_accel *acc = gsl_interp_accel_alloc();

gsl_spline_init(Y, &xvals[0], &yvals[0], yvals.size());

double y[N];

for (int i = 0; i < N; ++i) {
    double x = xmin + i*dx; // Where xmin is the smallest x value and dx
                            // is (xmax-xmin)/N
    y[i] = gsl_spline_eval(Y, x, acc);
}

然后在计算幂律和数据之间差异的函数中，我确保使用相同的 xmin 和 dx，以便符号中的 Yi 函数的x值相同。

struct data {
    size_t n;
    double *y;
    double xmin;
    double xmax;
};

int powerLaw(const gsl_vector *x, void *dat, gsl_vector *f) {
    size_t n = ((data *) dat)->n;
    double *y = ((data *) dat)->y;
    double xmin = ((data *) dat)->xmin;
    double xmax = ((data *) dat)->xmax;
    double dx = (xmax-xmin)/double(n);

    double A = gsl_vector_get(x, 0);
    double alpha = gsl_vector_get(x, 1);

    for (int i = 0; i < n; ++i) {
        double xval = xmin + double(i)*dx;
        double Yi = A*pow(xval,alpha);
        gsl_vector_set(f, i, Yi - y[i]);
    }

    return GSL_SUCCESS;
}

之后，来自 gnuplot 和 GSL 的值非常吻合，gnuplot 给出的幅度为 123.196 +/- 0.04484，指数为 -1.13275 +/- 0.001903，GSL 给出的值为 123.20464 +/- 0.98008 和 -1.13272 +/- 0.00707。拟合的结果如下图所示，其中 Fit 来自 gnuplot，g(x) 来自 GSL（注意：我不希望幂律与数据精确匹配，但对于我的目的）。gnuplot 和 GSL 的拟合几乎相同。

绘制来自 gnuplot 和 GSL 的数据和拟合图。

我会在对您的问题的评论中提到这一点，但是由于我从来没有在这里提出问题并且从未回答过任何问题，所以我没有足够的代表。

于 2016-07-13T18:08:16.730 回答