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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <mach.h>
include "curfitdef.h"
define COV Memr[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.
procedure cverrors (cv, y, w, yfit, npts, chisqr, errors)
pointer cv # curve descriptor
real y[ARB] # data points
real yfit[ARB] # fitted data points
real w[ARB] # array of weights
int npts # number of points
real chisqr # reduced chi-squared of fit
real errors[ARB] # errors in coefficients
int i, n, nfree
real variance, chisq, hold
pointer sp, covptr
begin
# allocate space for covariance vector
call smark (sp)
call salloc (covptr, CV_NCOEFF(cv), TY_REAL)
# 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 aclrr (COV(covptr), CV_NCOEFF(cv))
COV(covptr+i-1) = 1.
call rcvchoslv (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
|
