;+ ; NAME: ; REGRESS2 ; ; PURPOSE: ; Multiple linear regression fit. ; Fit the function: ; Y(i) = A(0)*X(0,i) + A(1)*X(1,i) + ... + ; A(Nterms-1)*X(Nterms-1,i) ; ; CATEGORY: ; G2 - Correlation and regression analysis. ; ; CALLING SEQUENCE: ; Result = REGRESS(X, Y, W [, YFIT, SIGMA, FTEST, R, RMUL, CHISQ]) ; ; INPUTS: ; X: array of independent variable data. X must ; be dimensioned (Nterms, Npoints) where there are Nterms ; coefficients to be found (independent variables) and ; Npoints of samples. ; ; Y: vector of dependent variable points, must have Npoints ; elements. ; ; W: vector of weights for each equation, must be a Npoints ; elements vector. For instrumental weighting ; w(i) = 1/standard_deviation(Y(i)), for statistical ; weighting w(i) = 1./Y(i). For no weighting set w(i)=1, ; and also set the RELATIVE_WEIGHT keyword. ; ; OUTPUTS: ; Function result = coefficients = vector of ; Nterms elements. Returned as a column vector. ; ; OPTIONAL OUTPUT PARAMETERS: ; Yfit: array of calculated values of Y, Npoints elements. ; ; Sigma: Vector of standard deviations for coefficients. ; ; Ftest: value of F for test of fit. ; ; Rmul: multiple linear correlation coefficient. ; ; R: Vector of linear correlation coefficient. ; ; Chisq: Reduced chi squared. ; ; KEYWORDS: ;RELATIVE_WEIGHT: if this keyword is non-zero, the input weights ; (W vector) are assumed to be relative values, and not based ; on known uncertainties in the Y vector. This is the case for ; no weighting W(*) = 1. ; ; PROCEDURE: ; Adapted from the program REGRES, Page 172, Bevington, Data Reduction ; and Error Analysis for the Physical Sciences, 1969. ; ; MODIFICATION HISTORY: ; Written, DMS, RSI, September, 1982. ; Added RELATIVE_WEIGHT keyword, W. Landsman, August 1991 ; 29-AUG-1994: H.C. Wen - Used simpler, clearer algorithm to determine ; fit coefficients. The constant term, A0 is now just one ; of the X(iterms,*) vectors, enabling the algorithm to ; now return the sigma associated with this constant term. ; I also made a special provision for the case when ; Nterm = 1; namely, "inverting" a 1x1 matrix, i.e. scalar. ; 26-MAR-1996: Added the DOUBLE and CHECK keywords to the call to DETERM. ; 02-APR-1996: Test matrix singularity using second argument in INVERT ; instead of call to DETERM. ;- function REGRESS2,X,Y,W,Yfit,Sigma,Ftest,R,Rmul,Chisq, RELATIVE_WEIGHT=relative_weight On_error,2 ;Return to caller if an error occurs NP = N_PARAMS() if (NP lt 3) or (NP gt 9) then \$ message,'Must be called with 3-9 parameters: '+\$ 'X, Y, W [, Yfit, Sigma, Ftest, R, RMul, Chisq]' ; Determine the length of these arrays and the number of sources SX = SIZE( X ) SY = SIZE( Y ) nterm = SX(1) npts = SY(1) if (N_ELEMENTS(W) NE SY(1)) OR \$ (SX(0) NE 2) OR (SY(1) NE SX(2)) THEN \$ message, 'Incompatible arrays.' WW = REPLICATE(1.,nterm) # W curv = ( X*WW ) # TRANSPOSE( X ) beta = X # (Y*W) if nterm eq 1 then begin sigma = 1./sqrt(curv) X_coeff= beta/curv endif else begin err = INVERT( curv, status ) if (status eq 1) then begin print,'det( Curvature matrix )=0 .. Using REGRESS' X1 = X linechk = X(0,0) - X(0,fix( npts*randomu(seed) )) if linechk eq 0 then begin print,'Cannot determine sigma for CONSTANT' X1 = X1(1:nterm-1,*) endif coeff = REGRESS( X1,Y,W,Yfit,A0, Sigma,Ftest,R,Rmul,Chisq) if linechk eq 0 then begin coeff = [A0,reform(coeff)] Sigma = [ 0,reform(Sigma)] R = [ 0,R] endif return, coeff endif else if (status eq 2) then begin print,'WARNING -- small pivot element used in matrix inversion.' print,' significant accuracy probably lost.' endif diag = indgen( nterm ) sigma = sqrt( err( diag,diag ) ) X_coeff = err # beta endelse Yfit = TRANSPOSE(X_coeff # X) dof = npts - nterm > 1 Chisq = TOTAL( (Y-Yfit)^2.*W ) Chisq = Chisq/dof ; To calculate the "test of fit" parameters, we revert back to the original ; cryptic routine in REGRESS1. Because the constant term (if any) is now ; included in the X variable, NPAR = NTERM_regress2 = NTERM_regress1 + 1. if nterm eq 1 then goto, SKIP SW = TOTAL(W) ;SUM OF WEIGHTS YMEAN = TOTAL(Y*W)/SW ;Y MEAN XMEAN = (X * (REPLICATE(1.,NTERM) # W)) # REPLICATE(1./SW,NPTS) WMEAN = SW/NPTS WW = W/WMEAN ; NFREE = NPTS-1 ;DEGS OF FREEDOM SIGMAY = SQRT(TOTAL(WW * (Y-YMEAN)^2)/NFREE) ;W*(Y(I)-YMEAN) XX = X- XMEAN # REPLICATE(1.,NPTS) ;X(J,I) - XMEAN(I) WX = REPLICATE(1.,NTERM) # WW * XX ;W(I)*(X(J,I)-XMEAN(I)) SIGMAX = SQRT( XX*WX # REPLICATE(1./NFREE,NPTS)) ;W(I)*(X(J,I)-XM)*(X(K,I)-XM) R = WX #(Y - YMEAN) / (SIGMAX * SIGMAY * NFREE) WW1 = WX # TRANSPOSE(XX) ARRAY = INVERT(WW1/(NFREE * SIGMAX #SIGMAX)) A = (R # ARRAY)*(SIGMAY/SIGMAX) ;GET COEFFICIENTS FREEN = NPTS-NTERM > 1 ;DEGS OF FREEDOM, AT LEAST 1. CHISQ = TOTAL(WW*(Y-YFIT)^2)*WMEAN/FREEN ;WEIGHTED CHI SQUARED IF KEYWORD_SET(relative_weight) then varnce = chisq \$ else varnce = 1./wmean RMUL = TOTAL(A*R*SIGMAX/SIGMAY) ;MULTIPLE LIN REG COEFF IF RMUL LT 1. THEN FTEST = RMUL/(NTERM-1)/ ((1.-RMUL)/FREEN) ELSE FTEST=1.E6 RMUL = SQRT(RMUL) SKIP: return, X_coeff end