; ; Computes variance versus n for densities ; FOR ADJUSTED TIMESERIES ; loadct,39 sc=findgen(8)*30.+40. & sc=sc*(!d.n_colors-1)/255. ss=[1,4,6,1,4,6,1,4] yr=[0.25,0.25,0.4,0.4,0.8,0.5,0.5,1.5] multi_plot,nrow=4 if !d.name eq 'X' then begin initx window,ysize=850 endif ; ; Get details of regions & rbar ; restore,filename='reglists.idlsave' restore,filename='rbar_dens_full.idlsave' ; ; Do each region separately ; for i = 0 , nreg-1 do begin ; ; Read in timeseries of densities ; restore,filename='densadj_'+regname(i)+'.idlsave' keeplist=where(x lt 1991,nkeep) y=densadj(keeplist) x=x(keeplist) ; ; Compute mean n (chrons), ncore (cores) and variance from overlapping 50 yrs ; nyr=n_elements(x) segsize=51 nout=nyr-segsize+1 runx=fltarr(nout) runvar=fltarr(nout) runn=fltarr(nout) runncore=fltarr(nout) for j = 0 , nyr-segsize do begin yseg=y(j:j+segsize-1) xseg=x(j:j+segsize-1) nseg=n(j:j+segsize-1) ncoreseg=neff(j:j+segsize-1) runx(j)=total(xseg)/total(finite(xseg)) runn(j)=total(nseg)/total(finite(yseg)) ; only mean n where we have data runncore(j)=total(ncoreseg)/total(finite(yseg)) keeplist=where(finite(yseg),nkeep) if nkeep lt segsize then begin runvar(j)=!values.f_nan endif else begin ygot=yseg(keeplist) dummy=moment(ygot) runvar(j)=dummy(1) endelse endfor ; ; Now plot the results ; maxn=max(runn,/nan,k) maxvar=runvar(k) pause plot,x,y,xrange=[1600,2000],yrange=[-3,3],$ title=regname(i),ytitle='max density timeseries' plot,runx,runvar,xrange=[1600,2000],yrange=[0.,yr(i)],$ ytitle='running variance' plot,runx,runn,xrange=[1600,2000],ytitle='no. of chronologies' plot,runn,runvar,xrange=[0,100],psym=1,yrange=[0,3],$ ytitle='running variance',xtitle='no. of chronologies' ; if i eq 0 then plot,(1.+(runn-1.)*allrbar(i))/runn,runvar,$ ; xrange=[0,0.6],psym=1,yrange=[0,6], $ ; xtitle='1 / effective mean no. of chronologies',/xstyle, $ ; ytitle='running variance' $ ; else oplot,(1.+(runn-1.)*allrbar(i))/runn,runvar,psym=ss(i),color=sc(i) ; Keep the data for n<= 50 for fitting a line if i eq 0 then begin ; For all region, n > 50 always. Just initialise the arrays fitx=[0.] fity=[0.] endif else begin keeplist=where(finite(runvar) and (runn le 50),nkeep) if nkeep gt 0 then begin fitx=[fitx,runn(keeplist)] fity=[fity,1./runvar(keeplist)] endif endelse endfor ;plot,[0,1],/nodata,xstyle=5,ystyle=5 ;legend,regname,textcolors=sc,colors=sc,psym=ss ; end