svdfit.cpp

c经典算法(包括源码)

void svdfit(double x[], double y[], double sig[], int ndata, double a[],
			int ma, double u[], double v[], double w[], int mp, int np,
			double& chisq, char* funcs)
{
	int i, j;
    double thresh,tmp,tol = 0.00001;
    double b[1001], afunc[51];
    for (i = 1; i<=ndata; i++)
	{
        if (!strcmp(funcs,"fpoly"))
		{
			fpoly(x[i], afunc, ma);
		}
        if (!strcmp(funcs,"fleg"))
		{
			fleg(x[i], afunc, ma);
		}
        tmp = 1.0 / sig[i];
        for (j = 1; j<=ma; j++)
		{
            u[(i-1)*ma+j] = afunc[j] * tmp;
        }
        b[i] = y[i] * tmp;
    }
    svdcmp(u, ndata, ma, w, v);
    double wmax = 0.0;
    for (j = 1; j<=ma; j++)
	{
        if (w[j] > wmax)
		{
			wmax = w[j];
		}
    }
    thresh = tol * wmax;
    for (j = 1; j<=ma; j++)
	{
        if (w[j] < thresh)
		{
			w[j] = 0.0;
		}
    }
    svbksb(u, w, v, ndata, ma, b, a);
    chisq = 0.0;
    for (i = 1; i<=ndata; i++)
	{
        if (!strcmp(funcs,"fpoly"))
		{
			fpoly(x[i], afunc, ma);
		}
        if (!strcmp(funcs,"fleg"))
		{
			fleg(x[i], afunc, ma);
		}
        double sum1 = 0.0;
        for (j = 1; j<=ma; j++)
		{
            sum1 = sum1 + a[j] * afunc[j];
        }
        chisq = chisq + ((y[i] - sum1) / sig[i]) * ((y[i] - sum1) / sig[i]);
    }
}