1、一元线性模型:
公式:Y=A+ B*X
IDL提供了两个函数来实现
LinFit:The LINFIT function fits the paired data {xi, yi} to the
linear model, y = A+ Bx, by minimizing the chi-square error
statistic.
函数如下:
FUNCTION LinFit, x, y, $
chisqr = chisqr, Double = DoubleIn, prob = prob,
$
sig_ab = sig_ab, sigma = sigma, $
COVAR=covar,YFIT=yfit,
MEASURE_ERRORS=measure_errors, $
SDEV = sdevIn ; obsolete
keyword (still works)
COMPILE_OPT idl2
ON_ERROR, 2
nX = N_ELEMENTS(x)
nY = N_ELEMENTS(y)
if nX ne nY then $
MESSAGE, "X
and Y must be vectors of equal length."
;If the DOUBLE keyword is not set then
the internal precision and
;result are identical to the type of
input.
double = (N_ELEMENTS(doubleIn) gt 0) ?
KEYWORD_SET(doubleIn) : $
(SIZE(x,
/TYPE) eq 5) || (SIZE(y, /TYPE) eq 5)
isSdev = N_ELEMENTS(sdevIn) GT 0
isMeasure = N_ELEMENTS(measure_errors) GT
0
IF (isSdev OR isMeasure) THEN BEGIN
IF (isSdev
AND isMeasure) THEN $
MESSAGE,'Conflicting keywords SDEV and MEASURE_ERRORS.'
sdev =
isMeasure ? measure_errors : sdevIn
nsdev =
N_ELEMENTS(sdev)
IF (nsdev NE
nX) THEN MESSAGE, $
'MEASURE_ERRORS must have the number of elements as X and Y.'
ENDIF ELSE BEGIN
sdev =
double ? 1d : 1.0
nsdev =
0
ENDELSE
; for explanation of constants see Numerical Recipes sec.
15-2
if nsdev eq nX then begin ;Standard deviations
are supplied.
wt = 1/(sdev^(double ? 2d : 2.0))
ss = TOTAL(wt, DOUBLE=double)
sx = TOTAL(wt * x, DOUBLE=double)
sy = TOTAL(wt * y, DOUBLE=double)
t
= (x - sx/ss) / sdev
b = TOTAL(t
* y / sdev, DOUBLE=double)
endif else begin
ss = nX
sx = TOTAL(x, DOUBLE=double)
sy = TOTAL(y, DOUBLE=double)
t = x -
sx/ss
b = TOTAL(t
* y)
endelse
st2 = TOTAL(t^2, DOUBLE=double)
IF (NOT double) THEN BEGIN
ss = FLOAT(ss)
sx = FLOAT(sx)
sy = FLOAT(sy)
st2 = FLOAT(st2)
b = FLOAT(b)
ENDIF
; parameter estimates
b = b / st2
a = (sy - sx * b) / ss
; error estimates
sdeva = SQRT((1.0 + sx * sx / (ss * st2)) /
ss)
sdevb = SQRT(1.0 / st2)
covar = -sx/(ss*st2)
covar = [[sdeva^2, covar], [covar,
sdevb^2]]
yfit = b*x + a
if nsdev ne 0 then begin
chisqr =
TOTAL( ((y - yfit) / sdev)^2, Double = Double )
if ~Double
then $
chisqr = FLOAT(chisqr)
prob = 1 -
IGAMMA(0.5*(nX-2), 0.5*chisqr)
endif else begin
chisqr =
TOTAL( (y - yfit)^2, Double = Double )
if ~Double
then $
chisqr = FLOAT(chisqr)
prob =
double ? 1d : 1.0
sdevdat =
SQRT(chisqr / (nX-2))
sdeva =
sdeva * sdevdat
sdevb =
sdevb * sdevdat
endelse
sigma = (sig_ab = [sdeva,
sdevb])
RETURN, [a, b]
END
LadFit:The LADFIT function fits the paired data {xi, yi} to the
linear model, y = A+ Bx, using a "robust" least absolute deviation
method.
函数如下:
FUNCTION MDfunc, b, x, y, a, absdev,
EPS=eps
COMPILE_OPT hidden
a = MEDIAN(y - b*x,
/EVEN) ;EVEN
has no effect with odd counts.
d = y - (b * x + a)
absdev = total(abs(d))
nz = where(y ne 0.0, nzcount)
if nzcount ne 0 then d[nz] = d[nz] / abs(y[nz])
;Normalize
nz = where(abs(d) gt EPS, nzcount)
if nzcount ne 0 then
$
;Sign fcn, +1 for d > 0, -1 for d <
0, else 0.
return,
total(x[nz] * ((d[nz] gt 0) - fix(d[nz] lt 0))) $
else return, 0.0
END
;-------------------------------------------------------------------------
FUNCTION LadFit, xIn, yIn, absdev = absdev, Double = DoubleIn
ON_ERROR, 2
nX = N_ELEMENTS(xIn)
if nX ne N_ELEMENTS(yIn) then $
MESSAGE, "X
and Y must be vectors of equal length."
;If the DOUBLE keyword is not set then
the internal precision and
;result are identical to the type of
input.
Double = (N_ELEMENTS(doubleIn) gt 0) ?
KEYWORD_SET(doubleIn) : $
(SIZE(xIn,/TYPE) eq 5 or SIZE(yIn,/TYPE) eq 5)
x = Double ? DOUBLE(xIn) :
FLOAT(xIn)
y = Double ? DOUBLE(yIn) : FLOAT(yIn)
sx = TOTAL(x, Double = Double)
sy = TOTAL(y, Double = Double)
; the variance computation is sensitive to
roundoff, so we do this
; math in DP
sxy = TOTAL(DOUBLE(x)*DOUBLE(y), Double =
Double)
sxx = TOTAL(DOUBLE(x)*DOUBLE(x), Double =
Double)
del = DOUBLE(nX) * sxx - sx^2
if (del eq 0.0) then
begin
;All X's are the same
result =
[MEDIAN(y, /EVEN), 0.0] ;Bisect the range w/ a flat line
return,
Double ? result : FLOAT(result)
endif
aa = (sxx * sy - sx * sxy) / del ;Least
squares solution y = x * aa + bb
bb = (nX * sxy - sx * sy) / del
chisqr = TOTAL((y - (aa + bb*x))^2, Double =
Double)
sigb = sqrt(chisqr /
del)
;Standard deviation
b1 = bb
eps = Double ? 1d-7 : 1e-7
f1 = MDfunc(b1, x, y, aa, absdev,
EPS=eps)
; Quick return. The initial least squares
gradient is the LAD solution.
if (f1 eq 0.) then begin
bb=b1
goto,
done
endif
delb = ((f1 ge 0) ? 3.0 : -3.0) *
sigb
b2 = b1 + delb
f2 = MDfunc(b2, x, y, aa, absdev,
EPS=eps)
while (f1*f2 gt 0) do
begin
;Bracket the zero of the function
b1 = b2
f1 = f2
b2 = b1 + delb
f2 = MDfunc(b2, x, y, aa, absdev, EPS=eps)
endwhile
; In case we finish early.
bb = b2
f = f2
;Narrow tolerance to refine 0 of
fcn.
sigb = (Double ? 0.01d : 0.01) * sigb
while (abs(b2-b1) gt sigb
&& f ne 0) do begin ;bisection of
interval b1,b2.
bb = 0.5 *
(b1 + b2)
if (bb eq b1
|| bb eq b2) then $
break
f =
MDfunc(bb, x, y, aa, absdev, EPS=eps)
if (f*f1 ge
0) then begin
f1 = f
b1 = bb
endif else
begin
f2 = f
b2 = bb
endelse
endwhile
done: absdev = absdev /
nX
RETURN, Double ? [aa, bb] : FLOAT([aa,
bb])
END
效果:LadFit的结果优于LinFit,但是使用LadFit首先需要对X进行升序排列
2、多元线性模型:
多元线性拟合公式:Y = A + B*X1+ C*X2 ......
Regress:The REGRESS function performs a multiple linear
regression fit and returns an Nterm-element column vector of
coefficients.
必须注意X1,X2,X3不能完全相同
3、一元多次模型:
公式:Y = A + B*X + C*X^2 + D*X^3 +.....
Poly_Fit:The POLY_FIT function performs a least-square
polynomial fit with optional weighting and returns a vector of
coefficients.
下面是作的一个线性拟合功能的界面,可以通过采样数据进行回归系数的计算,模型暂时设计了上述三种:
http://rsgisman.bokee.com/photo/view.fcgi?mode=3&id=8139318中的线性拟合模型" TITLE="[转载]IDL 中的线性拟合模型" />
另外IDL还提供了其他非线性拟合模型的函数,也许以后会用到,此次没有用到就不列举了。