1 files changed, 111 insertions, 0 deletions
diff --git a/math/nlfit/nlerrors.gx b/math/nlfit/nlerrors.gx
new file mode 100644
index 00000000..061d2214
--- /dev/null
+++ b/math/nlfit/nlerrors.gx
@@ -0,0 +1,111 @@
+include <mach.h>
+include <math/nlfit.h>
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+define	COV	Mem$t[P2P($1)]	# element of COV
+
+
+# NLERRORS -- Procedure to calculate the reduced chi-squared of the fit
+# and the errors 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.
+
+procedure nlerrors$t (nl, z, zfit, w, npts, variance, chisqr, errors)
+
+pointer	nl		# curve descriptor
+PIXEL	z[ARB]		# data points
+PIXEL	zfit[ARB]	# fitted data points
+PIXEL	w[ARB]		# array of weights
+int	npts		# number of points
+PIXEL	variance	# variance of the fit
+PIXEL	chisqr		# reduced chi-squared of fit (output)
+PIXEL	errors[ARB]	# errors in coefficients (output)
+
+int	i, n, nfree
+pointer	sp, covptr
+PIXEL   factor
+
+begin
+	# Allocate space for covariance vector.
+	call smark (sp)
+	call salloc (covptr, NL_NPARAMS(nl), TY_PIXEL)
+
+	# Estimate the variance and reduce chi-squared of the fit.
+	n = 0
+	variance = PIXEL (0.0)
+	chisqr = PIXEL (0.0)
+	do i = 1, npts {
+	    if (w[i] <= PIXEL (0.0))
+	        next
+	    factor = (z[i] - zfit[i]) ** 2
+	    variance = variance + factor
+	    chisqr = chisqr + factor * w[i]
+	    n = n + 1
+	}
+
+	# Calculate the reduced chi-squared.
+	nfree = n - NL_NFPARAMS(nl)
+	if (nfree  > 0) {
+	    variance = variance / nfree
+	    chisqr = chisqr / nfree
+	} else {
+	    variance = PIXEL (0.0)
+	    chisqr = PIXEL (0.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 and not the reduced chi-squared
+
+	if (abs (chisqr - variance) <= EPSILON$T) {
+	    if (nfree > 0)
+		factor = chisqr
+	    else
+		factor = PIXEL (0.0)
+	} else
+	    factor = 1.
+
+	# Calculate the  errors in the coefficients.
+	call aclr$t (errors, NL_NPARAMS(nl))
+	call nlinvert$t (ALPHA(NL_ALPHA(nl)), CHOFAC(NL_CHOFAC(nl)),
+	    COV(covptr), errors, PLIST(NL_PLIST(nl)), NL_NFPARAMS(nl), factor)
+
+	call sfree (sp)
+end
+
+
+# NLINVERT -- Procedure to invert matrix and compute errors
+
+procedure nlinvert$t (alpha, chofac, cov, errors, list, nfit, variance)
+
+PIXEL	alpha[nfit,ARB]		# data matrix
+PIXEL	chofac[nfit, ARB]	# cholesky factorization
+PIXEL	cov[ARB]		# covariance vector
+PIXEL	errors[ARB]		# computed errors
+int	list[ARB]		# list of active parameters
+int	nfit			# number of fitted parameters
+PIXEL	variance		# variance of the fit
+
+int	i, ier
+
+begin
+	# Factorize the data matrix to determine the errors.
+	call nl_chfac$t (alpha, nfit, nfit, chofac, ier)
+
+	# Estimate the errors.
+	do i = 1, nfit {
+	    call aclr$t (cov, nfit)
+	    cov[i] = PIXEL (1.0)
+	    call nl_chslv$t (chofac, nfit, nfit, cov, cov)
+	    if (cov[i] >= PIXEL (0.0))
+	        errors[list[i]] = sqrt (cov[i] * variance)
+	    else
+		errors[list[i]] = PIXEL (0.0)
+	}
+end