; 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