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还提供了其他非线性拟合模型的函数,也许以后会用到,此次没有用到就不列举了。