加载中…
个人资料
  • 博客等级:
  • 博客积分:
  • 博客访问:
  • 关注人气:
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
正文 字体大小:

[转载]IDL 中的线性拟合模型

(2014-12-04 11:01:53)
标签:

转载

分类: 遥感
原文地址:IDL 中的线性拟合模型作者:jpg_PNG

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

0

  

新浪BLOG意见反馈留言板 欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

新浪公司 版权所有