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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
|
c subroutine regres.f
c
c source
c Bevington, pages 172-175.
c
c purpose
c make a mulitple linear regression fit to data with a specified
c function which is linear in coefficients
c
c usage
c call regres (x, y, sigmay, npts, nterms, m, mode, yfit,
c a0, a, sigma0, sigmaa, r, rmul, chisqr, ftest)
c
c description of parameters
c x - array of points for independent variable
c y - array of points for dependent variable
c sigmay - array of standard deviations for y data points
c npts - number of pairs of data points
c nterms - number of coefficients
c m - array of inclusion/rejection criteria for fctn
c mode - determines method of weighting least-squares fit
c +1 (instrumental) weight(i) = 1./sigmay(i)**2
c 0 (no weighting) weight(i) = 1.
c -1 (statistical) weight(i) = 1./y(i)
c yfit - array of calculated values of y
c a0 - constant term
c a - array of coefficients
c sigma0 - standard deviation of a0
c sigmaa - array of standard deviations for coefficients
c r - array of linear correlation coefficients
c rmul - multiple linear correlation coefficient
c chisqr - reduced chi square for fit
c ftest - value of f for test of fit
c
c subroutines and function subprograms required
c fctn (x, i, j, m)
c evaluates the function for the jth term and the ith data point
c using the array m to specify terms in the function
c matinv (array, nterms, det)
c inverts a symmetric two-dimensional matrix of degree nterms
c and calculates its determinant
c
c comments
c (dim npts changed 100->1000 21-may-84 dct)
c dimension statement valid for npts up to 100 and nterms up to 10
c sigmaag changed to sigmaa in statement following statement 132
c
subroutine pfregs (x,y,sigmay,npts,nterms,m,mode,yfit,
*a0,a,sigma0,sigmaa,r,rmul,chisqr,ftest,fctn)
double precision array,sum,ymean,sigma,chisq,xmean,sigmax
dimension x(1),y(1),sigmay(1),m(1),yfit(1),a(1),sigmaa(1),
*r(1)
dimension weight(1000),xmean(10),sigmax(10),array(10,10)
external fctn
c
c initialize sums and arrays
c
11 sum=0.
ymean=0.
sigma=0.
chisq=0.
rmul=0.
do 17 i=1,npts
17 yfit(i)=0.
21 do 28 j=1,nterms
xmean(j)=0.
sigmax(j)=0.
r(j)=0.
a(j)=0.
sigmaa(j)=0.
do 28 k=1,nterms
28 array(j,k)=0.
c
c accumulate weighted sums
c
30 do 50 i=1,npts
31 if (mode) 32,37,39
32 if (y(i)) 35,37,33
33 weight(i)=1./y(i)
goto 41
35 weight(i)=1./(-y(i))
goto 41
37 weight(i)=1.
goto 41
39 weight(i)=1./sigmay(i)**2
41 sum=sum+weight(i)
ymean=ymean+weight(i)*y(i)
do 44 j=1,nterms
44 xmean(j)=xmean(j)+weight(i)*fctn(x,i,j,m)
50 continue
51 ymean=ymean/sum
do 53 j=1,nterms
53 xmean(j)=xmean(j)/sum
fnpts=npts
wmean=sum/fnpts
do 57 i=1,npts
57 weight(i)=weight(i)/wmean
c
c accumulate matrices r and array
c
61 do 67 i=1,npts
sigma=sigma+weight(i)*(y(i)-ymean)**2
do 67 j=1,nterms
sigmax(j)=sigmax(j)+weight(i)*(fctn(x,i,j,m)-xmean(j))**2
r(j)=r(j)+weight(i)*(fctn(x,i,j,m)-xmean(j))*(y(i)-ymean)
do 67 k=1,j
67 array(j,k)=array(j,k)+weight(i)*(fctn(x,i,j,m)-xmean(j))*
*(fctn(x,i,k,m)-xmean(k))
71 free1=npts-1
72 sigma=dsqrt(sigma/free1)
do 78 j=1,nterms
74 sigmax(j)=dsqrt(sigmax(j)/free1)
r(j)=r(j)/(free1*sigmax(j)*sigma)
do 78 k=1,j
array(j,k)=array(j,k)/(free1*sigmax(j)*sigmax(k))
78 array(k,j)=array(j,k)
c
c invert symmetric matrix
c
81 call matinv (array,nterms,det)
if (det) 101,91,101
91 a0=0.
sigma0=0.
rmul=0.
chisqr=0.
ftest=0.
goto 150
c
c calculate coefficients, fit, and chi square
c
101 a0=ymean
102 do 108 j=1,nterms
do 104 k=1,nterms
104 a(j)=a(j)+r(k)*array(j,k)
105 a(j)=a(j)*sigma/sigmax(j)
106 a0=a0-a(j)*xmean(j)
107 do 108 i=1,npts
108 yfit(i)=yfit(i)+a(j)*fctn(x,i,j,m)
111 do 113 i=1,npts
yfit(i)=yfit(i)+a0
113 chisq=chisq+weight(i)*(y(i)-yfit(i))**2
freen=npts-nterms-1
115 chisqr=chisq*wmean/freen
c
c calculate uncertainties
c
121 if (mode) 122,124,122
122 varnce=1./wmean
goto 131
124 varnce=chisqr
131 do 133 j=1,nterms
132 sigmaa(j)=array(j,j)*varnce/(free1*sigmax(j)**2)
if (sigmaa(j)) 835, 835, 836
835 sigmaa(j) = 0.0
goto 133
836 sigmaa(j)=sqrt(sigmaa(j))
133 rmul=rmul+a(j)*r(j)*sigmax(j)/sigma
freej=nterms
c +noao: When rmul = 1, the following division (stmt 135) would blow up.
c It has been changed so ftest is set to -99999. in this case.
if (rmul) 935, 136, 136
935 ftest = -99999.
rmul = -99999.
goto 141
c -noao
136 if (1.0 - abs(rmul)) 1035, 1036, 1037
1035 rmul=-99999.
ftest = -99999.
goto 141
1036 ftest = -99999.
rmul = 1.0
goto 141
1037 ftest=(rmul/freej)/((1.-rmul)/freen)
rmul=sqrt(rmul)
141 sigma0=varnce/fnpts
do 145 j=1,nterms
do 145 k=1,nterms
145 sigma0=sigma0+varnce*xmean(j)*xmean(k)*array(j,k)/
*(free1*sigmax(j)*sigmax(k))
146 sigma0=sqrt(sigma0)
150 return
end
|
