0

我正在尝试使用优化库 MPFIT 将高斯函数拟合到我的数据中。实际上，此代码是 MPFIT 库附带的示例代码的一部分。原始代码自动在内部以数字方式计算函数的导数，并且运行良好。MPFIT 库还允许用户提供函数衍生物。这就是问题开始的地方。这是用于计算残差和函数的一阶偏导数的函数。

int gaussfunc(int m, int n, double *p, double *dy, double **derivs, void *vars)
{
  int i,j;
  struct vars_struct *v = (struct vars_struct *) vars;
  double *x, *y, *ey;
  double a = p[1];
  double b = p[2];
  double c = p[3];
  double d = p[0];
  x = v->x;
  y = v->y;
  ey = v->ey;

  for (i=0; i<m; i++)
  {
    dy[i] = (y[i] - (a*exp(-(x[i]-b)*(x[i]-b)/(2*c*c))+d))/ey[i];
  }

  // the code below this point is the code I added to calculate the derivatives.
  if(derivs)
  {
      for(j = 0; j < n; j++)
      {
          if (derivs[j])
          {
              for (i = 0; i < m; i++)
              {
                  double da =     exp(-(x[i]-b)*(x[i]-b)/(2*c*c));
                  double db = a * exp(-(x[i]-b)*(x[i]-b)/(2*c*c)) * (x[i]-b)/(c*c);
                  double dc = a * exp(-(x[i]-b)*(x[i]-b)/(2*c*c)) * (x[i]-b)*(x[i]-b)/(c*c*c);
                  double dd = 1;
                  double foo;
                  if     (j == 0) foo = dd;
                  else if(j == 1) foo = da;
                  else if(j == 2) foo = db;
                  else if(j == 3) foo = dc;
                  derivs[j][i] = foo;
              }
          }
      }
  }

  return 0;
}

'if (derivs)' 行上方的代码是原始代码，但经过重构，下面的代码是我用于计算导数的代码。我相信我的数学是正确的，它们已通过https://math.stackexchange.com/questions/716545/calculating-the-first-order-partial-derivatives-of-the-gaussian-function/716553验证

有没有人在将 MPFIT 与用户定义的导数一起使用时遇到同样的问题？

谢谢你。

4

1 回答 1

0

因为残差的计算是(DATA-MODEL)/SIGMA，所以导数应该是：

[-d(模型)/d(参数)]/sigma.

所以，这一行：

derivs[j][i] = foo;

变成：

derivs[j][i] = -foo/ey[i];

问题解决了！谢谢！

于 2014-03-19T02:03:30.483 回答