!*************************stat5*********************************** ! statistical analysis using B Bunk 1993 ! o data blocking rev 6/2014 ! o autocorrelation analysis ! o jackknife ! nvar : no. of variables ! nbmax : max. no. of blocks per variable ! - - - - - - - - - - - - - - - - - - - - - - - - - - module stat5 public clear5, accum5, aver5, sigma5, var5, covar5, tau5, tauint5 public jackout5, jackeval5, jack5 public save5, savef5, get5, getf5 private integer, save :: nvar = 0, nbmax = 0 integer, save, dimension(:), allocatable :: nbl, lbl, lnew real(8), save, dimension(:,:), allocatable:: blksum real(8), save, dimension(:), allocatable :: sqsum contains ! - - - - - - - - - - - - - - - - - - - - - - - - - - ! allocate and clear counters subroutine clear5(nvar1,nbmax1) real(8) zero parameter(zero=0.) if(nvar1.lt.1) stop 'error in clear5: nvar is invalid' if(nbmax1.lt.2 .or. mod(nbmax1,2).ne.0) & stop 'error in clear5: nbmax is invalid' if (nvar1 /= nvar .or. nbmax1 /= nbmax) then if (nvar > 0) then deallocate(nbl, lbl, lnew, blksum, sqsum, stat=ierror) if(ierror /= 0) stop 'error in clear5: deallocate failed' endif allocate(nbl(nvar1), lbl(nvar1), lnew(nvar1), & blksum(nvar1,nbmax1), sqsum(nvar1), stat=ierror) if (ierror /= 0) stop 'error in clear5: allocate failed' endif nvar=nvar1 nbmax=nbmax1 do 10 ivar=1,nvar nbl(ivar)=0 lbl(ivar)=1 lnew(ivar)=0 sqsum(ivar)=zero do 11 ibl=1,nbmax blksum(ivar,ibl)=zero 11 continue 10 continue end subroutine clear5 ! - - - - - - - - - - - - - - - - - - - - - - - - - - ! accumulate data subroutine accum5(ivar,value) real(8) value,zero parameter(zero=0.) if(ivar.lt.1 .or. ivar.gt.nvar) & stop 'error in accum5: ivar out of range' if(nbl(ivar) .eq. nbmax)then do 10 ibl=1,nbmax/2 blksum(ivar,ibl)=blksum(ivar,2*ibl-1) + blksum(ivar,2*ibl) 10 continue do 20 ibl=nbmax/2+1,nbmax blksum(ivar,ibl)=zero 20 continue nbl(ivar)=nbmax/2 lbl(ivar)=2*lbl(ivar) endif iblnew=nbl(ivar) + 1 blksum(ivar,iblnew)=blksum(ivar,iblnew) + value sqsum(ivar)=sqsum(ivar) + value**2 lnew(ivar)=lnew(ivar) + 1 if(lnew(ivar) .eq. lbl(ivar))then nbl(ivar)=iblnew lnew(ivar)=0 endif end subroutine accum5 ! - - - - - - - - - - - - - - - - - - - - - - - - - - ! compute averages real(8) function aver5(ivar) real(8) zero parameter(zero=0.) if(ivar.lt.1 .or. ivar.gt.nvar) & stop 'error in aver5: ivar out of range' aver5=zero nmeas=nbl(ivar)*lbl(ivar) + lnew(ivar) if(nmeas.eq.0) return do 10 ibl=1,min(nbl(ivar)+1,nbmax) aver5=aver5 + blksum(ivar,ibl) 10 continue aver5=aver5/nmeas end function aver5 ! - - - - - - - - - - - - - - - - - - - - - - - - - - ! compute errors real(8) function sigma5(ivar) real(8) zero parameter(zero=0.) if(ivar.lt.1 .or. ivar.gt.nvar) & stop 'error in sigma5: ivar out of range' sigma5=sqrt(max( var5(ivar) , zero )) end function sigma5 ! - - - - - - - - - - - - - - - - - - - - - - - - - - ! compute variances real(8) function var5(ivar) real(8) av,d,gam,gamvar,zero,one dimension d(nbmax) parameter(zero=0.,one=1.) if(ivar.lt.1 .or. ivar.gt.nvar) & stop 'error in var5: ivar out of range' var5=zero nb=nbl(ivar) if(nb.lt.2)return nmeas=nb*lbl(ivar) + lnew(ivar) av=zero do 10 ib=1,nb av=av + blksum(ivar,ib) 10 continue av=av/nb do 20 ib=1,nb d(ib)=blksum(ivar,ib) - av var5=var5 + d(ib)**2 20 continue if(var5 .le. zero)return var5=var5/nb gamvar=var5**2 do 30 it=1,nb-1 gam=zero do 40 ib=1,nb-it gam=gam + d(ib)*d(ib+it) 40 continue gam=gam/(nb-it) if( gam .le. sqrt(gamvar/(nb-it)) )goto 31 var5=var5 + 2*gam gamvar=gamvar + 2*gam**2 30 continue 31 continue var5=(var5/nmeas)/lbl(ivar) end function var5 ! - - - - - - - - - - - - - - - - - - - - - - - - - - ! compute elements of the covariance matrix ! - sum off-diagonal gammas ! - can violate positivity ! better use covar5 ! real(8) function cov5(ivar,jvar) real(8) avi,avj,di,dj,gamij,gamji,gamii,gamjj,gamvar real(8) sgn,zero,one dimension di(nbmax),dj(nbmax) parameter(zero=0.,one=1.) if(ivar.lt.1 .or. ivar.gt.nvar) & stop 'error in cov5: ivar out of range' if(jvar.lt.1 .or. jvar.gt.nvar) & stop 'error in cov5: jvar out of range' if(ivar .eq. jvar)then cov5=var5(ivar) return endif cov5=zero nb=nbl(ivar) if(nb.lt.2)return nmeas=nb*lbl(ivar) + lnew(ivar) if(nmeas .ne. nbl(jvar)*lbl(jvar)+lnew(jvar) )then print*,'warning in cov5: mismatch of measurements' return endif avi=zero avj=zero do 10 ib=1,nb avi=avi + blksum(ivar,ib) avj=avj + blksum(jvar,ib) 10 continue avi=avi/nb avj=avj/nb gamii=zero gamjj=zero do 20 ib=1,nb di(ib)=blksum(ivar,ib) - avi dj(ib)=blksum(jvar,ib) - avj cov5=cov5 + di(ib)*dj(ib) gamii=gamii + di(ib)**2 gamjj=gamjj + dj(ib)**2 20 continue cov5=cov5/nb sgn=sign(one,cov5) gamvar=gamii*gamjj/nb**2 + cov5**2 do 30 it=1,nb-1 gamij=zero gamji=zero gamii=zero gamjj=zero do 40 ib=1,nb-it gamij=gamij + di(ib+it)*dj(ib) gamji=gamji + dj(ib+it)*di(ib) gamii=gamii + di(ib)*di(ib+it) gamjj=gamjj + dj(ib)*dj(ib+it) 40 continue gamij=gamij/(nb-it) gamji=gamji/(nb-it) gamii=gamii/(nb-it) gamjj=gamjj/(nb-it) if( sgn*(gamij+gamji) .le. sqrt(2*gamvar/(nb-it)) )goto 31 cov5=cov5 + gamij + gamji gamvar=gamvar + 2*(gamii*gamjj + gamij*gamji) 30 continue 31 continue cov5=(cov5/nmeas)/lbl(ivar) end function cov5 ! - - - - - - - - - - - - - - - - - - - - - - - - - - ! compute elements of the covariance matrix ! - rescale non-diagonal elements ! with autocorrelation factors ! from the diagonal ! - positive matrix real(8) function covar5(ivar,jvar) real(8) avi,avj,di,dj,vari,varj,zero parameter(zero=0.) if(ivar.lt.1 .or. ivar.gt.nvar) & stop 'error in covar5: ivar out of range' if(jvar.lt.1 .or. jvar.gt.nvar) & stop 'error in covar5: jvar out of range' if(ivar .eq. jvar)then covar5=var5(ivar) return endif covar5=zero nb=nbl(ivar) if(nb.lt.2)return nmeas=nb*lbl(ivar) + lnew(ivar) if(nmeas .ne. nbl(jvar)*lbl(jvar)+lnew(jvar) )then print*,'warning in covar5: mismatch of measurements' return endif avi=zero avj=zero do 10 ib=1,nb avi=avi + blksum(ivar,ib) avj=avj + blksum(jvar,ib) 10 continue avi=avi/nb avj=avj/nb vari=zero varj=zero do 20 ib=1,nb di=blksum(ivar,ib) - avi dj=blksum(jvar,ib) - avj covar5=covar5 + di*dj vari=vari + di**2 varj=varj + dj**2 20 continue if(vari.le.zero .or. varj.le.zero)return covar5=covar5 * sqrt( var5(ivar)*var5(jvar)/(vari*varj) ) end function covar5 ! - - - - - - - - - - - - - - - - - - - - - - - - - - ! estimate integrated auto-correlations: ! var(average) = var(single)/nmeas * rsq ! rsq = coth( 1/(2*tau) ) ! = 2 * tauint real(8) function tau5(ivar) real(8) var,rsq,zero,one,oneeps parameter(zero=0.,one=1.,oneeps=1.000001d0) if(ivar.lt.1 .or. ivar.gt.nvar) & stop 'error in tau5: ivar out of range' tau5=zero nmeas=nbl(ivar)*lbl(ivar) + lnew(ivar) if(nmeas.lt.2) return var=sqsum(ivar)/nmeas - aver5(ivar)**2 if(var.le.zero) return rsq=var5(ivar) / var * nmeas if(rsq.le.oneeps) return tau5=one/log( (rsq+one)/(rsq-one) ) end function tau5 ! . . . . . . . . . . . . . . . . . . . . . . . . . . real(8) function rsq5(ivar) real(8) var,rsq,zero,one parameter(zero=0.,one=1.) if(ivar.lt.1 .or. ivar.gt.nvar) & stop 'error in rsq5: ivar out of range' rsq5=one nmeas=nbl(ivar)*lbl(ivar) + lnew(ivar) if(nmeas.lt.2) return var=sqsum(ivar)/nmeas - aver5(ivar)**2 if(var.le.zero) return rsq5=var5(ivar) / var * nmeas end function rsq5 ! . . . . . . . . . . . . . . . . . . . . . . . . . . real(8) function tauint5(ivar) real(8) var,rsq,zero,one parameter(zero=0.,one=1.) if(ivar.lt.1 .or. ivar.gt.nvar) & stop 'error in tauint5: ivar out of range' tauint5=.5 nmeas=nbl(ivar)*lbl(ivar) + lnew(ivar) if(nmeas.lt.2) return var=sqsum(ivar)/nmeas - aver5(ivar)**2 if(var.le.zero) return tauint5=.5 * var5(ivar) / var * nmeas end function tauint5 ! - - - - - - - - - - - - - - - - - - - - - - - - - - ! compute jackknife blocks ! ivar : variable (input) ! nb : no. of blocks (output) ! bj : real(8) vector of jackknife blocks ! bj(ib), ib=1..nb (output) subroutine jackout5(ivar,nb,bj) real(8) bj dimension bj(nbmax) real(8) bsum,zero,one parameter(zero=0.,one=1.) if(ivar.lt.1 .or. ivar.gt.nvar) & stop 'error in jackout5: ivar out of range' nb=nbl(ivar) if(nb .lt. 2) stop 'error in jackout5: nb < 2' bsum=zero do 10 ib=1,nb bsum=bsum + blksum(ivar,ib) 10 continue do 20 ib=1,nb bj(ib)=(bsum - blksum(ivar,ib))/(lbl(ivar)*(nb-1)) 20 continue end subroutine jackout5 ! - - - - - - - - - - - - - - - - - - - - - - - - - - ! jackknife analysis for a vector of (function) values ! nb : number of jackknife values, nb > 1 (input) ! fj : real(8) vector of jackknife values, ! fj(ib), ib=1..nb (input) ! aver : function average (output) ! sigma : error estimate (output) subroutine jackeval5(nb,fj,aver,sigma) real(8) fj,aver,sigma dimension fj(nb) real(8) bsum,d,var,gam,gamvar,zero,one dimension d(nb) parameter(zero=0.,one=1.) aver=zero sigma=zero if(nb.lt.2)return bsum=zero do 10 ib=1,nb bsum=bsum + fj(ib) 10 continue aver=bsum/nb gam=zero do 20 ib=1,nb d(ib)=fj(ib) - aver gam=gam + d(ib)**2 20 continue if(gam .le. zero)return gam=gam/nb var=gam gamvar=gam**2 do 30 it=1,nb-1 gam=zero do 40 ib=1,nb-it gam=gam + d(ib)*d(ib+it) 40 continue gam=gam/(nb-it) if( gam .le. sqrt(gamvar/(nb-it)) )goto 31 var=var + 2*gam gamvar=gamvar + 2*gam**2 30 continue 31 continue sigma=(nb-1)*sqrt(var/nb) end subroutine jackeval5 ! - - - - - - - - - - - - - - - - - - - - - - - - - - ! jackknife evaluation - custom call ! fct : function of statistical variables (user-defined) ! function fct(nvar, a) ! real(8) a(nvar) ! with a(i), i=1..nvar : variables, as in stat5 ! aver : function average (output) ! sigma : error estimate (output) ! calls jackeval5 ! note: is is assumed that all (relevant) variables were accumulated ! in sync, with ivar=1 as prototype. subroutine jack5(fct,aver,sigma) real(8) fct,aver,sigma real(8) bsum,bj,fj,zero,one dimension bsum(nvar),bj(nvar),fj(nbmax) parameter(zero=0.,one=1.) external fct aver=zero sigma=zero nb=nbl(1) lb=lbl(1) if(nb.lt.2)return do 10 ivar=1,nvar bsum(ivar)=zero do 11 ib=1,nb bsum(ivar)=bsum(ivar) + blksum(ivar,ib) 11 continue 10 continue do 15 ib=1,nb do 16 ivar=1,nvar bj(ivar)=(bsum(ivar) - blksum(ivar,ib))/(lb*(nb-1)) 16 continue fj(ib)=fct(nvar,bj) 15 continue call jackeval5(nb,fj,aver,sigma) end subroutine jack5 ! - - - - - - - - - - - - - - - - - - - - - - - - - - ! save and restore counters (unformatted/formatted) subroutine save5(iunit) write(iunit)nvar,nbmax write(iunit)(nbl(i),lbl(i),lnew(i),i=1,nvar) & ,((blksum(i,ibl),ibl=1,nbmax),i=1,nvar) & ,(sqsum(i),i=1,nvar) end subroutine save5 subroutine savef5(iunit) write(iunit,*)nvar,nbmax write(iunit,*)(nbl(i),lbl(i),lnew(i),i=1,nvar) & ,((blksum(i,ibl),ibl=1,nbmax),i=1,nvar) & ,(sqsum(i),i=1,nvar) end subroutine savef5 subroutine get5(iunit) read(iunit)nvar1,nbmax1 if(nvar1.lt.1) stop 'error in get5: nvar is invalid' if(nbmax1.lt.2 .or. mod(nbmax1,2).ne.0) & stop 'error in get5: nbmax is invalid' if (nvar1 /= nvar .or. nbmax1 /= nbmax) then if (nvar > 0) then deallocate(nbl, lbl, lnew, blksum, sqsum, stat=ierror) if(ierror /= 0) stop 'error in get5: deallocate failed' endif allocate(nbl(nvar1), lbl(nvar1), lnew(nvar1), & blksum(nvar1,nbmax1), sqsum(nvar1), stat=ierror) if (ierror /= 0) stop 'error in get5: allocate failed' endif nvar=nvar1 nbmax=nbmax1 read(iunit)(nbl(i),lbl(i),lnew(i),i=1,nvar) & ,((blksum(i,ibl),ibl=1,nbmax),i=1,nvar) & ,(sqsum(i),i=1,nvar) end subroutine get5 subroutine getf5(iunit) read(iunit,*)nvar1,nbmax1 if(nvar1.lt.1) stop 'error in getf5: nvar is invalid' if(nbmax1.lt.2 .or. mod(nbmax1,2).ne.0) & stop 'error in getf5: nbmax is invalid' if (nvar1 /= nvar .or. nbmax1 /= nbmax) then if (nvar > 0) then deallocate(nbl, lbl, lnew, blksum, sqsum, stat=ierror) if(ierror /= 0) stop 'error in getf5: deallocate failed' endif allocate(nbl(nvar1), lbl(nvar1), lnew(nvar1), & blksum(nvar1,nbmax1), sqsum(nvar1), stat=ierror) if (ierror /= 0) stop 'error in getf5: allocate failed' endif nvar=nvar1 nbmax=nbmax1 read(iunit,*)(nbl(i),lbl(i),lnew(i),i=1,nvar) & ,((blksum(i,ibl),ibl=1,nbmax),i=1,nvar) & ,(sqsum(i),i=1,nvar) end subroutine getf5 end module stat5