2

我尝试根据数值配方使用 Levenberg-marquardt 的方法来拟合模型。问题是：它不会收敛，或者当它收敛时，它并不精确......或者至少协变矩阵很奇怪。

        int i=0;
        for (i = 0; i < 3e4; i++) {
            mrqmin(x, y, sig, NPCalib, a, ia, 3, covar, alpha, &chisk, afunc,
                &alamda);
            if (chisk < 1e-8)
                sumchisk++;
            if (sumchisk > 5)
                break;
            if (alamda > 1e8)
                alamda = 1e8;
        }

(x,y) 是 3 个点（双精度），可以很好地与 . 形式配合使用y=a(x-x0)^2。
像这样使用 sumchisk 是使用此功能的数值配方的建议。
alamda 在这里被封顶，否则可能会溢出。

其他定义和数据点：

double a[4] = {0.0, 0.0001, 100.0, -1};
int ia[4] = {0.0, 1, 1, 0};
double *x = {0.0, 799.157549545577, 799.92196995454, 800.683769692575};
double *y = {0.0, 524.26491, 525.26768, 526.26586};
double *sig = {0.0, 0.1*y[1], 0.1*y[2], 0.1*y[3]};
double **covar = new double*[4];    
covar[1] = new double[4];
covar[2] = new double[4];
covar[3] = new double[4];
double **alpha = new double*[4];    
alpha[1] = new double[4];
alpha[2] = new double[4];
alpha[3] = new double[4];
double chisk = 0;
double alamda = -1;

void afunc(int i, double x[], double a[], double *y, double dyda[], int ma)
{
    *y = a[1] * pow(x[i] + a[2], 2) / pow(1 + a[3] * CT[i - 1], 2);

    dyda[1] = pow(x[i] + a[2], 2) / pow(1 + a[3] * CT[i - 1], 2);
    dyda[2] = (2 * a[1] * (x[i] + a[2])) / pow
        (1 + a[3] * CalibTurn[i - 1], 2);
    dyda[3] = (-2 * a[1] * CT[i - 1] * pow(x[i] + a[2], 2)) / pow
        (1 + a[3] * CT[i - 1], 3);
}

我将 nr-sourcecode 更改为使用 double 而不是 float。没有使用第一个数组元素，因为它来自 fortran 代码，我不想更改这么小的细节。

该模型还包含一个 3. 参数，该参数未在此拟合中使用，因此保持为 a[3]=-1，因为 ia[3]=0。ia[]=1 表示参数即将拟合...

但是，现在我遇到的问题是有时这不会收敛。它以alamda=1e8和结束i=3e4。尤其是当我将门槛设置为chisk更低时。
参数集似乎很好，虽然...... chisk 大约是 1e-6 并且参数看起来很好，但是查看协变矩阵的对角线（应该给出每个参数的平方标准偏差），对于参数 0.0001，有一些像 ~800000 这样的垃圾。

有谁知道我在使用这个算法时做错了什么？开始时我需要写入 covar/alpha 的任何具体内容吗？sig可以这样设置吗？

4

0 回答 0