Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.stsci.edu/~robberto/Main/Software/IDL4pipeline/linearize_mr.pro Дата изменения: Thu Feb 27 21:56:49 2014 Дата индексирования: Sun Mar 2 09:21:03 2014 Кодировка: Поисковые слова: m 2

; docformat = 'IDL'
;+
; GLS linear fitter
;-
function linearize_mr,time,ramp
;this function implements the Generalized Least Square Fit to a signal ramp
;
;make the ramp parameters available to the public
common ramp_params,N,M,Ms,G,tf,tg,RON,rate
common graphics, DISPLAY
;
;define the y variables, referred to the the first read
delta_y = ramp - ramp[0];
; a0=linfit(time,ramp)
; d_y=ramp-(a0[0]+a0[1]*time)
;
;set the COVARIANCE MATRIX
SIGMA=DBLarr(N,N)
;
;to estimate the variance I need an initial estimate of b, either from the minmax...
b = Delta_y[N-1] / ( (N-1)*tg)
;or from the OLS fit
;..a0=linfit(time,ramp)
;b=a0[1]
;
;I calculate the first row/col of the covariance matrix out of the loop, just for convenience
SIGMA[0,0] = b*tf*(M+1)/(2.*M) + b*tf*1/(3.*M)*(M-1)*(M-2) + RON^2/M + (G*M/SQRT(12))^2 & $
SIGMA[0,1:*] = (M+1)/2.*b*tf
SIGMA[1:*,0] = (M+1)/2.*b*tf
;
;the other terms of the matrix should be under control...
for i=1,N-1 do begin & $
; b = Delta_y[i] / ((i-1)*tg + (m+1)/2.*tf)
; SIGMA[i,i] = b*(i-1)*tg + (2*M^2-3*M+7)/(6*M)*b*tf + RON^2/M + (G*M/SQRT(12))^2 & $ ;Robberto Eq. 16
SIGMA[i,i] = b*i*tg + b*tf*(M+1)/(2.*M) + b*tf*1/(3.*M)*(M-1)*(M-2) + RON^2/M + (G*M/SQRT(12))^2 & $ ;Robberto Eq. 16
IF I LT N-1 THEN BEGIN & $;diagonal terms
SIGMA[I,I+1:N-1] = REPLICATE(b*i*tg+b*(m+1)/2.*tf,N-(i+1)) & $ ;Robberto Eq.19
SIGMA[I+1:N-1,I] = REPLICATE(b*i*tg+b*(m+1)/2.*tf,N-(i+1)) & $
ENDIF & $
endfor
;stop
;
;INVERSE COVARIANCE (AKA WEIGHTS if DIAGONAL)
W_Y_1 = INVERT(SIGMA)
;
;A MATRIX
AA=fltarr(2,N)
AA[0,*]=1.
AA[1,*]=time
;
;W_A MATRIX
W_A=INVERT(TRANSPOSE(AA)##W_Y_1##AA)
;
;and here is we have the RESULTS:
eta=W_A##TRANSPOSE(AA)##W_Y_1##Delta_Y
;..eta=W_A##TRANSPOSE(AA)##W_Y_1##D_Y
;
;IN READABLE FORM
intercept = eta[0] + ramp[0]
;intercept = eta[0] + a0[0]
slope = eta[1]; no a0[1] term since I have not done a subtraction of a0[1]*time
;slope = a0[1]+eta[1]
;
;and here are the errors
delta_intercept=SQRT(W_A[0,0])
delta_slope=SQRT(W_A[1,1])
;
;if DISPLAY has been set, we check individual ellipses
IF N_ELEMENTS(DISPLAY) NE 0 THEN BEGIN
params=covar2ellipse(W_A)
plot,indgen(2),xrange=[-0.5,0.5]+rate,yrange=[-50,50]+10000.,/nodata,XTITLE='Slope', YTITLE='Intercept'
oplot,[rate],[1E4],psym=2,symsize=1.5
params=covar2ellipse(W_A)
oplot,[slope],[intercept],psym=4
tvellipse, params.major, params.minor, slope, intercept, -params.ANGLE
;NOTE: I am playing a little trick here to have the ellipse properly oriented and sized;
;should be done scaling the covariance, but this is for illustration purposes only...
;the conventional OLS fit gives
a0=linfit(time,ramp,SIGMA=da0)
oploterror,[a0[1]],[a0[0]],[da0[1]],[da0[0]],hatlen=3
stop
;
ENDIF
;
;plot ellipse
;=======================================================================
params=covar2ellipse(W_A)
plot,indgen(2),xrange=[slope-delta_slope*3,slope+delta_slope*3],yrange=[intercept-delta_intercept*3,intercept+delta_intercept*3],/nodata,XTITLE='Slope', YTITLE='Intercept'
tvellipse, params.major, params.minor, slope, intercept, -params.ANGLE, Color='yellow', /FILL
oplot,[rate],[1E4],psym=2,symsize=1.5,Color='black';,symsize=0.5
params=covar2ellipse(W_A)
oplot,[slope],[intercept],psym=4,Color='black',symsize=0.5
;the conventional OLS fit gives
a0=linfit(time,ramp,SIGMA=da0)
tvellipse,da0[1],da0[0], a0[1], a0[0], 0, Color='red', /FILL
oplot,[a0[1]],[a0[0]],color='black',psym=4,symsize=0.5
; stop
; oploterror,[a0[1]],[a0[0]],[da0[1]],[da0[0]],hatlen=3,color='black'
;=======================================================================
;
;RETURN THE ARRAY WITH THE FOUR RESULTS
return,[intercept,slope,delta_intercept,delta_slope]
end