pro mkcalibratemulti,por,pand,fitcoeff=fitcoeff,autocoeff=autocoeff,$ nlag=nlag ; ; Given a set of predictor time series (por[n,t]) and a predictand time series ; (pand[t]), with no missing data, multiple linear (least-squares) ; regression is performed to calibrate the predictors. ; fitcoeff is returned containing {r,alpha,beta,SEalpha,SEbeta,sig,rms}, ; where r=correlation(prediction,pand), alpha and beta are the offset and ; scaling coefficients to be applied to the predictors, their standard errors ; are also given, the significance (as a % confidence) is returned (taking ; into account autocorrelation), and also the root-mean-squared error after ; calibration. autocoeff is returned containing {lag1r of por, lag1r of pand, ; and lag1r of residuals, effective independent sample size}. ; ; When computing the signficance of the correlation, the number of lagged ; autocorrelations to use can be specified in nlag (default is 1). This ; is usually only needed when using smoothed series that no longer look ; like simple red noise (when nlag=9 typically works well). ; The lag1r of each predictor is computed, and then only the maximum of these ; is used (to reduce sample size), although all are returned. ; ; NB. The multiple regression function does not compute a standard error ; for the constant (intercept) coefficient, so this is set to missing. ; ; Check time series lengths and number of predictors ; psize=size(por) case psize(0) of 1: begin nval=psize(1) npred=1 por=reform(por,npred,nval) end 2: begin nval=psize(2) npred=psize(1) end else: message,'por[n,t] must be 1D or 2D' endcase if n_elements(pand) ne nval then message,'Incompatible time series' if n_elements(nlag) eq 0 then nlag=1 ;print,npred,nval ; ; Compute multiple regression ; wts=fltarr(nval)+1. bcoeff=regress(por,pand,wts,/relative_weight,predval,acoeff,bse,fvalue,$ sepr,mulr) bcoeff=reform(bcoeff) ;print,bcoeff ;print,sepr ;print,fvalue ; ; Compute residuals resids=pand-predval ; ; Compute RMS error sqerr=resids^2 rmse=sqrt( total(sqerr)/nval ) ; ; Compute lag-1 autocorrelations alag=indgen(nlag)+1 ac1all=fltarr(npred,nlag) for ipred = 0 , npred-1 do begin ac1all(ipred,*)=a_correlate(reform(por(ipred,*)),alag) endfor dummy=max(reform(ac1all(*,0)),maxele) ;print,dummy ac1por=reform(ac1all(maxele,*)) ac1pand=a_correlate(pand,alag) ac1resid=a_correlate(resids,alag) ; ; Compute significance of r(por,pand), taking autocorrelation into account rr = (ac1por(0) > 0.) * (ac1pand(0) > 0.) if nlag eq 1 then ess=(nval-2.)*(1.-rr)/(1.+rr)+2. $ else ess=(nval-2.)/(1.+2.*total(rr))+2. v1=npred v2=(ess-npred-1) > 1 tsig=f_pdf(fvalue,v1,v2) tsig=tsig*100. ; fitcoeff = { r: mulr, alpha: acoeff, beta: bcoeff, sealpha: !values.f_nan,$ sebeta: bse, sig: tsig, rms: rmse } autocoeff= { ac1por: reform(ac1all(*,0)), ac1pand: ac1pand(0), $ ac1resid: ac1resid(0), ess:ess } ; return end