/usr/share/gnudatalanguage/astrolib/sixlin.pro is in gdl-astrolib 2018.02.16+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 | pro sixlin,xx,yy,a,siga,b,sigb,weight=weight
;+
; NAME:
; SIXLIN
; PURPOSE:
; Compute linear regression coefficients by six different methods.
; EXPLANATION:
; Adapted from the FORTRAN program (Rev. 1.1) supplied by Isobe,
; Feigelson, Akritas, and Babu Ap. J. Vol. 364, p. 104 (1990).
; Suggested when there is no understanding about the nature of the
; scatter about a linear relation, and NOT when the errors in the
; variable are calculable.
;
; CALLING SEQUENCE:
; SIXLIN, xx, yy, a, siga, b, sigb, [WEIGHT = ]
;
; INPUTS:
; XX - vector of X values
; YY - vector of Y values, same number of elements as XX
;
; OUTPUTS:
; A - Vector of 6 Y intercept coefficients
; SIGA - Vector of standard deviations of 6 Y intercepts
; B - Vector of 6 slope coefficients
; SIGB - Vector of standard deviations of slope coefficients
;
; The output variables are computed using linear regression for each of
; the following 6 cases:
; (0) Ordinary Least Squares (OLS) Y vs. X (c.f. linfit.pro)
; (1) Ordinary Least Squares X vs. Y
; (2) Ordinary Least Squares Bisector
; (3) Orthogonal Reduced Major Axis
; (4) Reduced Major-Axis
; (5) Mean ordinary Least Squares
;
; OPTIONAL INPUT KEYWORD:
; WEIGHT - vector of weights, same number of elements as XX and YY
; For 1 sigma Gausssian errors, the weights are 1/sigma^2 but
; the weight vector can be more general. Default is no
; weighting.
; NOTES:
; Isobe et al. make the following recommendations
;
; (1) If the different linear regression methods yield similar results
; then quoting OLS(Y|X) is probably the most familiar.
;
; (2) If the linear relation is to be used to predict Y vs. X then
; OLS(Y|X) should be used.
;
; (3) If the goal is to determine the functional relationship between
; X and Y then the OLS bisector is recommended.
;
; REVISION HISTORY:
; Written Wayne Landsman February, 1991
; Corrected sigma calculations February, 1992
; Added WEIGHT keyword J. Moustakas February 2007
;-
compile_opt idl2
On_error, 2 ;Return to Caller
if N_params() LT 5 then begin
print,'Syntax - SIXLIN, xx, yy, a, siga, b, sigb, {WEIGHT =]'
return
endif
b = dblarr(6) & siga = b & sigb =b
x = double(xx) ;Keep input X and Y vectors unmodified
y = double(yy)
rn = N_elements(x)
if rn LT 2 then $
message,'Input X and Y vectors must contain at least 2 data points'
if rn NE N_elements(y) then $
message,'Input X and Y vectors must contain equal number of data points'
if (n_elements(weight) eq 0L) then weight = replicate(1.0,rn) else begin
if (rn ne n_elements(weight)) then $
message,'Input X and WEIGHT vectors must contain equal number of data points'
endelse
; Compute averages and sums
sumw = total(weight)
xavg = total( weight * x)/sumw
yavg = total( weight * y)/sumw
x = x - xavg
y = y - yavg
sxx = total( weight * x^2)
syy = total( weight * y^2)
sxy = total( weight * x*y)
if sxy EQ 0. then $
message,'SXY is zero, SIXLIN is terminated'
if sxy LT 0. then sign = -1.0 else sign = 1.0
; Compute the slope coefficients
b[0] = sxy / sxx
b[1] = syy / sxy
b[2] = (b[0]*b[1] - 1.D + sqrt((1.D + b[0]^2)*(1.D +b[1]^2)))/(b[0] + b[1] )
b[3] = 0.5 * ( b[1] - 1.D/b[0] + sign*sqrt(4.0D + (b[1]-1.0/b[0])^2))
b[4] = sign*sqrt( b[0]*b[1] )
b[5] = 0.5 * ( b[0] + b[1] )
; Compute Intercept Coefficients
a = yavg - b*xavg
; Prepare for computation of variances
gam1 = b[2] / ( (b[0] + b[1]) * $
sqrt( (1.D + b[0]^2)*(1.D + b[1]^2)) )
gam2 = b[3] / (sqrt( 4.D*b[0]^2 + ( b[0]*b[1] - 1.D)^2))
sum1 = total( weight * ( x*( y - b[0]*x ) )^2)
sum2 = total( weight * ( y*( y - b[1]*x ) )^2)
sum3 = total( weight * x * y * ( y - b[0]*x) * (y - b[1]*x ) )
cov = sum3 / ( b[0]*sxx^2 )
; Compute variances of the slope coefficients
sigb[0] = sum1 / sxx^2
sigb[1] = sum2 / sxy^2
sigb[2] = (gam1^2) * ( ( (1.D + b[1]^2) ^2 )*sigb[0] + $
2.D*(1.D + b[0]^2) * (1.D + b[1]^2)*cov + $
( (1.D + b[0]^2)^2)*sigb[1] )
sigb[3] = (gam2^2)*( sigb[0]/b[0]^2 + 2.D*cov + b[0]^2*sigb[1] )
sigb[4] = 0.25*(b[1]*sigb[1]/b[1] + $
2.D*cov + b[0]*sigb[1]/b[1] )
sigb[5] = 0.25*(sigb[0] + 2.D*cov + sigb[1] )
; Compute variances of the intercept coefficients
siga[0] = total( weight * ( ( y - b[0]*x) * (1.D - sumw*xavg*x/sxx) )^2 )
siga[1] = total( weight * ( ( y - b[1]*x) * (1.D - sumw*xavg*y/sxy) )^2 )
siga[2] = total( weight * ( (x * (y - b[0]*x) * (1.D + b[1]^2) / sxx + $
y * (y - b[1]*x) * (1.D + b[0]^2) / sxy)* $
gam1 * xavg * sumw - y + b[2] * x) ^ 2)
siga[3] = total( weight * ( ( x * ( y - b[0]*x) / sxx + $
y * ( y - b[1]*x) * b[0]^2/ sxy) * gam2 * $
xavg * sumw / sqrt( b[0]^2) - y + b[3]*x) ^ 2 )
siga[4] = total( weight * ( ( x * ( y - b[0] * x) * sqrt( b[1] / b[0] ) / sxx + $
y * ( y - b[1] * x) * sqrt( b[0] / b[1] ) / sxy) * $
0.5 * sumw * xavg - y + b[4] * x)^2 )
siga[5] = total( weight * ( (x * ( y - b[0] * x) / sxx + $
y * ( y - b[1] * x) / sxy)* $
0.5 * sumw * xavg - y + b[5]*x )^2 )
; Convert variances to standard deviation
sigb = sqrt(sigb)
siga = sqrt(siga)/sumw
return
end
|