mrqmin.cpp

提供了几种数据拟合的c++代码：1 直线拟合(FIT)；2 线性最小二乘法((LFIT, COVSRT), (SVDFIT, SVDVAR)(oddity), (FPOLY, FLEG)(exampl

void mrqmin(double x[], double y[], double sig[], int ndata, double a[],
			int ma, int lista[], int mfit, double covar[], double alpha[],
			int nca, double& chisq, double& alamda)
{
	static kk,i,j,k,ihit;
    static double atry[21], beta[21], da[21];
	static double ochisq;
    if (alamda < 0.0)
	{
        kk = mfit + 1;
        for (j = 1; j<=ma; j++)
		{
            ihit = 0;
            for (k = 1; k<=mfit; k++)
			{
                if (lista[k] == j)
				{
					ihit = ihit + 1;
				}
            }
            if (ihit == 0)
			{
                lista[kk] = j;
                kk = kk + 1;
			}
            else
			{
				if (ihit > 1)
				{
					cout<< "improper permutation in lista"<<endl;
					return;
				}
            }
        }
        if (kk != ma + 1)
		{
			cout<< "improper permutation in lista"<<endl;
		}
        alamda = 0.001;
        mrqcof(x, y, sig, ndata, a, ma, lista, mfit, alpha, beta, nca, chisq);
        ochisq = chisq;
        for (j = 1; j<=ma; j++)
		{
            atry[j] = a[j];
        }
    }
    for (j = 1; j<=mfit; j++)
	{
        for (k = 1; k<=mfit; k++)
		{
            covar[(j-1)*nca+k] = alpha[(j-1)*nca+k];
        }
        covar[(j-1)*nca+j] = alpha[(j-1)*nca+j] * (1.0 + alamda);
        da[j] = beta[j];
    }
    gaussj(covar, mfit, nca, da);
    if (alamda == 0.0)
	{
        covsrt(covar, nca, ma, lista, mfit);
		for (i=1; i<=20; i++)
		{
			da[i]=0;
			atry[i]=0;
		}
        return;
    }
    for (j = 1; j<=mfit; j++)
	{
        atry[lista[j]] = a[lista[j]] + da[j];
    }
    mrqcof(x, y, sig, ndata, atry, ma, lista, mfit, covar, da, nca, chisq);
    if (chisq < ochisq)
	{
        alamda = 0.1 * alamda;
        ochisq = chisq;
        for (j = 1; j<=mfit; j++)
		{
            for (k = 1; k<=mfit; k++)
			{
                alpha[(j-1)*nca+k] = covar[(j-1)*nca+k];
            }
            beta[j] = da[j];
            a[lista[j]] = atry[lista[j]];
        }
	}
    else
	{
        alamda = 10.0 * alamda;
        chisq = ochisq;
    }
	for (i=1; i<=20; i++)
	{
		da[i]=0;
		atry[i]=0;
	}
}