reg_multlin

hi,
I'm seeing some strange results when using reg_multlin to perform
multiple linear regression.

I'm calling reg_multlin multiple times (once at each grid point of a
73x144 grid). At most gridpoints the results look reasonable, but at a
handful of points I get results which are many orders of magnitude
larger (x 1.e20 or so). The strange thing is that when I add an unused
dummy array (test1 in the attachment, currently commented out) to my
script, the bad results go away and all looks good.

Has anyone had similar experiences using reg_multlin?
Thanks,
Gary

;=======================================================================
;
; Description: (1) Regress psi onto rmm1 and rmm2.
;
; Execution: ncl reg.ncl
;=======================================================================
;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;
begin

;***********************************************************
; read in psi data
;***********************************************************
 fileanom0 = "/Volumes/Disk/gbates/mjo/psi250.anom0.nc"
 file0 = addfile( fileanom0,"r")
 printVarSummary(file0)

 time = file0->time
 t_len = dimsizes(time)
 print(t_len)

 year = time / 1000000
 month = time / 10000 - year*100
 day = time/100 - year*10000 - month*100

 lat = file0->lat
 lat_len = dimsizes(lat)

 lon = file0->lon
 lon_len = dimsizes(lon)

 nl=12
 psi0 = file0->psi250
 printVarSummary(psi0)

 file2 = "/Volumes/Disk/gbates/mjo/psi250.anom"+nl+".nc"
 filef = addfile( file2,"r")
 printVarSummary(filef)

 psif = filef->psi250

;test1 = psi0

;***********************************************************
; read in normalized data
;***********************************************************
 fileorig = "/Volumes/Disk/gbates/mjo/3var.normanom.nc"
 file1 = addfile( fileorig,"r") ; open CCM

 var = file1->var
 printVarSummary(var)

;***********************************************************
; read in eofs
;***********************************************************
 fileeof = "/Volumes/Disk/gbates/mjo/3var.eof.nc"
 filee = addfile( fileeof,"r")

 eof=filee->eof
 printVarSummary(eof)

 eof1 = -1.*eof(0,0,:)
 eof2 = -1.*eof(0,1,:)

;***********************************************************
; Compute normalization factors and project eofs
; on data to get rmm1 and rmm2.
;***********************************************************
 var2 = var(lat|0,lon|:,time|:)
 varrun = runave( var(lat|0,lon|:,time|:), 120, 0 )
 var2 = var2 - varrun

 norm1 = eof1 # var2
 norm2 = eof2 # var2

 rmm1 = norm1 / dim_stddev( norm1 )
 rmm2 = norm2 / dim_stddev( norm2 )
 rmm1 = linmsg(rmm1,0)
 rmm2 = linmsg(rmm2,0)

;***********************************************************
; Set up needed arrays.
;***********************************************************
 rmm1bar = (/ -1.53, 0., 1.54, 0. /)
 rmm2bar = (/ 0., -1.69, 0., 1.61 /)

 p1 = new( (/lat_len, lon_len/), float)
 p1!0 = "lat"
 p1!1 = "lon"
 p1&lat = lat
 p1&lon = lon

 p2 = new( (/lat_len, lon_len/), float)
 p2!0 = "lat"
 p2!1 = "lon"
 p2&lat = lat
 p2&lon = lon

 x = new( (/3, t_len/), float)
 x@_FillValue= -999.9

 x(0,:) = 1.
 x(1,:) = rmm1(0:t_len-1)
 x(2,:) = rmm2(0:t_len-1)
 xstd = dim_stddev(x)
 print(xstd)

  
 delpsi = psi0
 printVarSummary(delpsi)
 
 np=0
 
   kount = 0
   delpsi = -999.9

   ndays=14
   do nt=0,t_len-ndays-1
     if(.not.ismissing(rmm1(nt)))
       if(month(nt).le.3 .or. month(nt).ge.11) then
         xabs=abs(rmm1(nt))
         yabs=abs(rmm2(nt))
         mag = sqrt(xabs*xabs+yabs*yabs)
         if(mag.gt.1.0) then
           delpsi(nt,:,:) = psif(nt,:,:) - psi0(nt+nl,:,:)
         end if
       end if
     end if
   end do

;***********************************************************
; Multiple linear regression.
;***********************************************************
 y = new( (/t_len/), float)
 y@_FillValue= -999.9

 do nlat=0,lat_len-1
 do nlon=0,lon_len-1
   y = delpsi(:,nlat,nlon) * 1.e-7

   beta = reg_multlin (y,x,False)

   b = beta*xstd ; standard regression coefficients

   p1(nlat,nlon) = (/b(1)/)
   p2(nlat,nlon) = (/b(2)/)
   if(b(1).gt.1000.) then
     print(b(1))
   end if
 end do
 end do

 p1 = rmm1bar(np)*p1 + rmm2bar(np)*p2
 printMinMax(p1, True)

end

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Dec 19 2007 - 10:26:31 MST