Initial commit

1 files changed, 91 insertions, 0 deletions

diff --git a/math/curfit/cverrors.gx b/math/curfit/cverrors.gx
new file mode 100644
index 00000000..07288c7f
--- /dev/null
+++ b/math/curfit/cverrors.gx
@@ -0,0 +1,91 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <mach.h>
+
+$if (datatype == r)
+include "curfitdef.h"
+$else
+include "dcurfitdef.h"
+$endif
+
+define	COV	Mem$t[P2P($1)]	# element of COV
+
+# CVERRORS -- Procedure to calculate the reduced chi-squared of the fit
+# and the standard deviations of the coefficients. First the variance
+# and the reduced chi-squared of the fit are estimated. If these two
+# quantities are identical the variance is used to scale the errors
+# in the coefficients. The errors in the coefficients are proportional
+# to the inverse diagonal elements of MATRIX.
+
+$if (datatype == r)
+procedure cverrors (cv, y, w, yfit, npts, chisqr, errors)
+$else
+procedure dcverrors (cv, y, w, yfit, npts, chisqr, errors)
+$endif
+
+pointer	cv		# curve descriptor
+PIXEL	y[ARB]		# data points
+PIXEL	yfit[ARB]	# fitted data points
+PIXEL	w[ARB]		# array of weights
+int	npts		# number of points
+PIXEL	chisqr		# reduced chi-squared of fit
+PIXEL	errors[ARB]	# errors in coefficients
+
+int	i, n, nfree
+PIXEL   variance, chisq, hold
+pointer	sp, covptr
+
+begin
+	# allocate space for covariance vector
+	call smark (sp)
+	call salloc (covptr, CV_NCOEFF(cv), TY_PIXEL)
+
+	# estimate the variance and chi-squared of the fit
+	n = 0
+	variance = 0.
+	chisq = 0.
+	do i = 1, npts {
+	    if (w[i] <= 0.0)
+		next
+	    hold = (y[i] - yfit[i]) ** 2
+	    variance = variance + hold
+	    chisq = chisq + hold * w[i]
+	    n = n + 1
+	}
+
+	# calculate the reduced chi-squared
+	nfree = n - CV_NCOEFF(cv)
+	if (nfree  > 0)
+	    chisqr = chisq / nfree
+	else
+	    chisqr = 0.
+
+	# if the variance equals the reduced chi_squared as in the
+	# case of uniform weights then scale the errors in the coefficients
+	# by the variance not the reduced chi-squared
+	if (abs (chisq - variance) <= DELTA)
+	    if (nfree > 0)
+	        variance = chisq / nfree
+	    else
+		variance = 0.
+	else
+	    variance = 1.
+
+	# calculate the  errors in the coefficients
+	# the inverse of MATRIX is calculated column by column
+	# the error of the j-th coefficient is the j-th element of the
+	# j-th column of the inverse matrix
+	do i = 1, CV_NCOEFF(cv) {
+	    call aclr$t (COV(covptr), CV_NCOEFF(cv))
+	    COV(covptr+i-1) = 1.
+	    call $tcvchoslv (CHOFAC(CV_CHOFAC(cv)), CV_ORDER(cv), CV_NCOEFF(cv),
+	    	    COV(covptr), COV(covptr))
+	    if (COV(covptr+i-1) >= 0.)
+	        errors[i] = sqrt (COV(covptr+i-1) * variance)
+	    else
+		errors[i] = 0.
+	}
+
+
+	call sfree (sp)
+end