readmann =function(k, url= "d:/climate/data/mann.2008/proxy/A_allproxyoriginal.zip",id=idorig) { handle <- unz(url,paste(idorig$id[k],"ppd",sep=".") ,"r") #formerly V4.M02 Data <- scan(handle);length(Data) Data=t(array(Data,dim=c(3,length(Data)/3))) Data=data.frame(Data[4:nrow(Data),]) names(Data)=c("year","proxy","count") close(handle) readmann=Data readmann} #strip #FUNCTION TO LISTFILES AND REMOVE SUFFIXES strip<-function(h,n) { h<-unlist(strsplit(h, "\\.")) h<-array(h,dim=c(n,length(h)/n)) h<-t(h) strip<-h[,1] strip } filter.combine.pad=function(x,a,M=NA,method="norm") { temp<-!is.na(x) w<-x[temp] N<-length(w) if(is.na(M)) M<-trunc (length(a)/2)+1 if(method=="norm") w<-c (rep(mean(w[1:M]),M),w,rep(mean(w[(N-M+1):N]),M) ) #if(method=="reflect") w<-c ( w[M:1],w,w[M:(N-M+1):N]) if(method=="mannomatic") {w<-c( c (w[1]-(w[(M+1):2]-w[1])),w, c(w[N]- (w[(N-1):(N-M)]-w[N]) ) ) } if(method=="mann2") w<-c (rep(w[1],M),w,rep(w[N],M) ) y<-filter(w,a,method="convolution",sides=2) y<-y[(M+1):(length(y)-M)] z<- rep(NA,length(x)) z[temp]<-y if(class(x)=="ts") z<-ts(z,start=tsp(x)[1],freq=tsp(x)[3]) filter.combine.pad<-cbind(x,z) #filter.combine<-y[,2] filter.combine.pad } ##gaussian.filter.weights gaussian.filter.weights<-function(N,year) { #N number of points; year - is number of years sigma<-year/6 i<-( (-(N-1)/2):((N-1)/2) ) /sigma w<- ((2*pi)^-0.5) * exp(-0.5* i^2) gaussian.filter.weights<-w/sum(w) gaussian.filter.weights } #GAUSSIAN FILTER #see http://www.ltrr.arizona.edu/~dmeko/notes_8.pdf #The appropriate Gaussian distribution therefore has standard deviation # s = ? / 6, # where ? is the desired wavelength at which the amplitude of frequency response is 0.5. The # filter weights are obtained by sampling the pdf of the standard normal distribution at # values i/s where i={-(N-1)/s, ..., 0, .. (N-1)/s } # MEKo has values to verify for case of N=21, year=10 #this is verified ##truncated.gauss.weights truncated.gauss.weights<-function(year) { a<-gaussian.filter.weights(2*year+1,year) temp<-(a>0.05*max(a)) a<-a[temp] a<-a/sum(a) truncated.gauss.weights<-a truncated.gauss.weights } ##JONES #JONES LOOKUP jones<-function(lat,long) { i<-18-floor(lat/5); j<- 36+ceiling(long/5); jones<-72*(i-1) + j jones } ##jones.inv jones.inv<-function(k) { j<- k%%72 temp<-k-j i<-(temp/72)+1 long<-(j-36)*5-2.5 lat<-(18-i)*5+2.5 jones.inv<-c(lat,long) jones.inv } #http://www.climateaudit.org/?p=4038#comment-304170 filtfiltx = function( b, a, x = NULL ) { if (is.null(x)) { x = a a = b$a b = b$b } len = length(x) nfilt = max(length(b),length(a)) nfact = 3*(nfilt-1) # padding size if (len <= nfact) stop("Data must have length more than 3 times filter order.") # Not enough data b = c(b, rep(0,nfilt))[1:nfilt] # zero-pad shorter of b and a a = c(a, rep(0,nfilt))[1:nfilt] # ditto y = c(2*x[1]-rev(x[2:(nfact+1)]), x, 2*x[len]-rev(x[(len-nfact):(len-1)])) y = rev(y[1] + filter(b,a,y-y[1])) y = rev(y[1] + filter(b,a,y-y[1])) y[(nfact+1):(len+nfact)] } mannsmooth=function(x,bwf=bf,M=5,method="filtfiltx") { #filtfiltx mannsmooth = x temp=!is.na(x);N=sum(temp) x=x[temp]; m0=mean(x) if (method=="center") x=x-m0 mannextend=array(rep(x,9),dim=c(N,9)) dummy1=cbind(rep(mean(x),M),x[M:1],2*x[1]-x[M:1]); dummy1=cbind(dummy1,dummy1,dummy1) #these are 3 padding cases at start: mean; y-reflection; xy-reflection (endpoint is center) dummy2=array(c(rep(mean(x), 3*M),rep(x[(N-1):(N-M)],3),rep(2*x[N]-x[(N-1):(N-M)],3)),dim=c(M,9)); #three padding cases at end mannextend=rbind(dummy1,mannextend,dummy2);dim(mannextend) #2023 9 mannextend=mannextend[,c(1,4,7,2,5,8,3,6,9)] #to match iconb, icone order working=apply(mannextend,2, function(x) filtfiltx(bwf,x) ) #this uses Roman M's implementation if(method=="filtfilt") working=apply(mannextend,2, function(x) filtfilt(bwf,x) ) test=working[(M+1):(N+M),] residssq=apply(test-x,2,function(x) sum(x^2)) #var in original but doesn't matter index=max((1:9)[residssq==min(residssq)]);index #the last value is chosen #tie breaker in Mann is the last series mannsmooth[temp]=test[,index] if(method=="center") mannsmooth=mannsmooth+m0 mannsmooth } #this extracts an underlying decadal series from a linear smoothed time series extractf=function(x) { r1=range(time(x)) year=c(time(x)) y=diff(x) z=round(diff(y),2) temp=!(z==0)&!is.na(z) time0=sort(c(r1,time(z)[temp]-1)) rx=round(x,2) temp=!is.na(match(year, time0));sum(temp) test=data.frame(time(x)[temp],rx[temp]/units) names(test)=c("year","proxy") #sum(!is.na(test$proxy[test$year>=1850])) #9 extractf=test# list(data=test,x=x,y=y,z=z,r1=r1) extractf } ##k=202;units=1;x=units*g(mannorig[[details$old[k] ]]) #collects according to X$year cumulative to end irrmatch=function(A,y) { gridded=data.frame(c(time(y)),y);names(gridded)=c("year","target") h= function( x) if(x[1]<=max(A$year)) min(A$year[ A$year>=x[1] ] ) else NA gridded$d= apply(gridded,1,h) working=tapply(gridded$target,gridded$d,mean,na.rm=T) test=match(A$year,as.numeric(names(working)));temp=!is.na(test) A$grid=NA;A$grid[temp]=working[test[temp]] irrmatch=A irrmatch} #rescales rescale=function (x,target, index=1850:1995) (x- mean(x[index-tsp(x)[1]+1],na.rm=T))* sd(target[index-tsp(target)[1]+1])/sd(x[index-tsp(x)[1]+1],na.rm=T)+mean(target[index-tsp(target)[1]+1]) #4. FUNCTION TO RETURN VERIFICATION PERIOD STATISTICS #calibration and verificaiton are both vectors of length 2 with start and end period #estimator and observed are time series which may have different start and end periods verification.stats<-function(estimator,observed=sparse,calibration=c(1850,1949),verification=c(1950,1995)) { combine<-ts.union(estimator,observed) index.cal<- (calibration[1]-tsp(combine)[1]+1):(calibration[2]-tsp(combine)[1]+1) index.ver<-(verification[1]-tsp(combine)[1]+1):(verification[2]-tsp(combine)[1]+1) mean.obs.cal<-mean ( combine[index.cal,"observed"],na.rm=TRUE ) #xmean mean.obs.ver<-mean( combine[index.ver,"observed"],na.rm=TRUE ) # mx cov.ver<-cov ( combine[index.ver,],use=use0) #was test.ver cov.cal<-cov(combine[index.cal,],use=use0) #was test R2.cal<-(cov.cal[1,2]*cov.cal[2,1])/(cov.cal[1,1]*cov.cal[2,2]) R2.ver<-(cov.ver[1,2]*cov.ver[2,1])/(cov.ver[1,1]*cov.ver[2,2]) RE.cal<-1-sum( (combine[index.cal,"estimator"]- combine[index.cal,"observed"]) ^2,na.rm=TRUE)/sum( (mean.obs.cal- combine[index.cal,"observed"]) ^2,na.rm=TRUE) RE.ver<-1-sum( (combine[index.ver,"estimator"]- combine[index.ver,"observed"]) ^2,na.rm=TRUE)/sum( (mean.obs.cal- combine[index.ver,"observed"]) ^2,na.rm=TRUE) CE<- 1-sum( (combine[index.ver,"estimator"]- combine[index.ver,"observed"]) ^2,na.rm=TRUE)/sum( (mean.obs.ver- combine[index.ver,"observed"]) ^2,na.rm=TRUE) verification.stats<-data.frame(RE.cal=RE.cal,RE.ver=RE.ver,R2.cal=R2.cal,R2.ver=R2.ver,CE=CE)# verification.stats } mannstats= function(estimator,observed=sparse,calibration=c(1850,1949),verification=c(1950,1995)) { combine<-ts.union(estimator,observed); index.cal<- (calibration[1]-tsp(combine)[1]+1):(calibration[2]-tsp(combine)[1]+1) index.ver<-(verification[1]-tsp(combine)[1]+1):(verification[2]-tsp(combine)[1]+1) combine=data.frame(combine) mean.obs.cal<-mean ( combine$observed[index.cal],na.rm=TRUE ) #xmean mean.obs.ver<-mean( combine$observed[index.ver],na.rm=TRUE ) # mx resid.ver= sum((combine$estimator[index.ver]-combine$observed[index.ver])^2) ssx.ver = sum( combine$observed[index.ver]^2) ############## mean is NOT subtracted for some reason betamulti=1-resid.ver/ssx.ver rcal=cor(combine$observed[index.cal],combine$estimator[index.cal],use="pairwise.complete.obs") rver=cor(combine$observed[index.ver],combine$estimator[index.ver],use="pairwise.complete.obs") R2.cal=rcal^2; R2.ver=rver^2 RE.ver=1-resid.ver/sum( (combine$estimator[index.ver]-mean.obs.cal)^2) CE= 1-resid.ver/sum( (combine$estimator[index.ver]-mean.obs.ver)^2) mannstats<-data.frame(RE.ver=RE.ver,R2.cal=R2.cal,R2.ver=R2.ver,CE=CE)# mannstats } #mannstats(x,target) betamulti=1