lfit.cpp

c经典算法(包括源码)

void lfit(double x[], double y[], double sig[], int ndata, double a[],int ma,
		  int lista[], int mfit, double covar[], int ncvm, double& chisq)
{
	int i,j,k,kk,ihit;
	double ym,sig2i,wt,sum;
    double beta[51], afunc[51];
    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 set in lista"<<endl;
			}
        }
    }
    if (kk != (ma + 1))
	{
		cout<<" improper set in lista"<<endl;
	}
    for (j = 1; j<=mfit; j++)
	{
        for (k = 1; k<=mfit; k++)
		{
            covar[(j-1)*ma+k] = 0.0;
        }
        beta[j] = 0.0;
    }
    for (i = 1; i<=ndata; i++)
	{
        funcs(x[i], afunc, ma);
        ym = y[i];
        if (mfit < ma)
		{
            for (j = mfit + 1; j<=ma; j++)
			{
                ym = ym - a[lista[j]] * afunc[lista[j]];
            }
        }
        sig2i = 1.0 / (sig[i] * sig[i]);
        for (j = 1; j<=mfit; j++)
		{
            wt = afunc[lista[j]] * sig2i;
            for (k = 1; k<=j; k++)
			{
                covar[(j-1)*ma+k] = covar[(j-1)*ma+k] + wt * afunc[lista[k]];
			}
            beta[j] = beta[j] + ym * wt;
        }
    }
    if (mfit > 1)
	{
        for (j = 2; j<=mfit; j++)
		{
            for (k = 1; k<=j - 1; k++)
			{
                covar[(k-1)*ma+j] = covar[(j-1)*ma+k];
            }
        }
    }
    gaussj(covar, mfit,ma,beta);
    for (j = 1; j<=mfit; j++)
	{
        a[lista[j]] = beta[j];
    }
    chisq = 0.0;
    for (i = 1; i<=ndata; i++)
	{
        funcs(x[i], afunc, ma);
        sum = 0.0;
        for (j = 1; j<=ma; j++)
		{
            sum = sum + a[j] * afunc[j];
        }
        chisq = chisq + ((y[i] - sum) / sig[i]) * ((y[i] - sum) / sig[i]);
    }
    covsrt(covar, ncvm, ma, lista, mfit);
}