From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001
From: Joseph Hunkeler <jhunkeler@gmail.com>
Date: Wed, 8 Jul 2015 20:46:52 -0400
Subject: Initial commit

---
 math/gsurfit/gserrors.gx | 90 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 90 insertions(+)
 create mode 100644 math/gsurfit/gserrors.gx

(limited to 'math/gsurfit/gserrors.gx')

diff --git a/math/gsurfit/gserrors.gx b/math/gsurfit/gserrors.gx
new file mode 100644
index 00000000..5f78cfab
--- /dev/null
+++ b/math/gsurfit/gserrors.gx
@@ -0,0 +1,90 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <mach.h>
+$if (datatype == r)
+include "gsurfitdef.h"
+$else
+include "dgsurfitdef.h"
+$endif
+
+define	COV	Mem$t[P2P($1)]	# element of COV
+
+# GSERRORS -- 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 gserrors (sf, z, w, zfit, chisqr, errors)
+$else
+procedure dgserrors (sf, z, w, zfit, chisqr, errors)
+$endif
+
+pointer	sf		# curve descriptor
+PIXEL	z[ARB]		# data points
+PIXEL	w[ARB]		# array of weights
+PIXEL	zfit[ARB]	# fitted data points
+PIXEL	chisqr		# reduced chi-squared of fit
+PIXEL	errors[ARB]	# errors in coefficients
+
+int	i, nfree
+PIXEL	variance, chisq, hold
+pointer	sp, covptr
+
+begin
+	# allocate space for covariance vector
+	call smark (sp)
+	$if (datatype == r)
+	call salloc (covptr, GS_NCOEFF(sf), TY_REAL)
+	$else
+	call salloc (covptr, GS_NCOEFF(sf), TY_DOUBLE)
+	$endif
+
+	# estimate the variance and chi-squared of the fit
+	variance = 0.
+	chisq = 0.
+	do i = 1, GS_NPTS(sf) {
+	    hold = (z[i] - zfit[i]) ** 2
+	    variance = variance + hold
+	    chisq = chisq + hold * w[i]
+	}
+
+	# calculate the reduced chi-squared
+	nfree = GS_NPTS(sf) - GS_NCOEFF(sf)
+	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, GS_NCOEFF(sf) {
+	    call aclr$t (COV(covptr), GS_NCOEFF(sf))
+	    COV(covptr+i-1) = 1.
+	    call $tgschoslv (CHOFAC(GS_CHOFAC(sf)), GS_NCOEFF(sf),
+	        GS_NCOEFF(sf), 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
-- 
cgit