1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
|
include <mach.h>
include <math/nlfit.h>
include "nlfitdefr.h"
define COV Memr[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 nlerrorsr (nl, z, zfit, w, npts, variance, chisqr, errors)
pointer nl # curve descriptor
real z[ARB] # data points
real zfit[ARB] # fitted data points
real w[ARB] # array of weights
int npts # number of points
real variance # variance of the fit
real chisqr # reduced chi-squared of fit (output)
real errors[ARB] # errors in coefficients (output)
int i, n, nfree
pointer sp, covptr
real factor
begin
# Allocate space for covariance vector.
call smark (sp)
call salloc (covptr, NL_NPARAMS(nl), TY_REAL)
# Estimate the variance and reduce chi-squared of the fit.
n = 0
variance = real (0.0)
chisqr = real (0.0)
do i = 1, npts {
if (w[i] <= real (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 = real (0.0)
chisqr = real (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) <= EPSILONR) {
if (nfree > 0)
factor = chisqr
else
factor = real (0.0)
} else
factor = 1.
# Calculate the errors in the coefficients.
call aclrr (errors, NL_NPARAMS(nl))
call nlinvertr (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 nlinvertr (alpha, chofac, cov, errors, list, nfit, variance)
real alpha[nfit,ARB] # data matrix
real chofac[nfit, ARB] # cholesky factorization
real cov[ARB] # covariance vector
real errors[ARB] # computed errors
int list[ARB] # list of active parameters
int nfit # number of fitted parameters
real variance # variance of the fit
int i, ier
begin
# Factorize the data matrix to determine the errors.
call nl_chfacr (alpha, nfit, nfit, chofac, ier)
# Estimate the errors.
do i = 1, nfit {
call aclrr (cov, nfit)
cov[i] = real (1.0)
call nl_chslvr (chofac, nfit, nfit, cov, cov)
if (cov[i] >= real (0.0))
errors[list[i]] = sqrt (cov[i] * variance)
else
errors[list[i]] = real (0.0)
}
end
|
