fit.cpp

此程序为VC++常用数值算法这本书附赠的光盘中包含了本书中全部的源代码

void fit(double x[], double y[], int ndata, double sig[], int mwt, double& a,
		double& b, double& siga, double& sigb, double& chi2, double& q)
{
	int i;
    double sigdat,t,sxoss,wt,ss,sx = 0.0;
    double sy = 0.0;
    double st2 = 0.0;
    b = 0.0;
    if (mwt != 0)
	{
        ss = 0.0;
        for (i = 1; i<=ndata; i++)
		{
            wt = 1.0 / (sig[i] * sig[i]);
            ss = ss + wt;
            sx = sx + x[i] * wt;
            sy = sy + y[i] * wt;
        }
	}
    else
	{
        for (i = 1; i<=ndata; i++)
		{
            sx = sx + x[i];
            sy = sy + y[i];
        }
        ss = ndata;
    }
    sxoss = sx / ss;
    if (mwt != 0)
	{
        for (i = 1; i<=ndata; i++)
		{
            t = (x[i] - sxoss) / sig[i];
            st2 = st2 + t * t;
            b = b + t * y[i] / sig[i];
        }
	}
    else
	{
        for (i = 1; i<=ndata; i++)
		{
            t = x[i] - sxoss;
            st2 = st2 + t * t;
            b = b + t * y[i];
		}
    }
    b = b / st2;
    a = (sy - sx * b) / ss;
    siga = sqrt((1.0 + sx * sx / (ss * st2)) / ss);
    sigb = sqrt(1.0 / st2);
    chi2 = 0.0;
    if (mwt == 0)
	{
        for (i = 1; i<=ndata; i++)
		{
            chi2 = chi2 + pow((y[i] - a - b * x[i]) , 2);
        }
        q = 1.0;
        sigdat = sqrt(chi2 / (ndata - 2));
        siga = siga * sigdat;
        sigb = sigb * sigdat;
	}
    else
	{
        for (i = 1; i<=ndata; i++)
		{
            chi2 = chi2 + pow(((y[i] - a - b * x[i]) / sig[i]) , 2);
        }
        q = gammq(0.5 * (ndata - 2), 0.5 * chi2);
    }
}