GPS
filenum=c(11:15,19:24)roads=c("鞍山西道","白堤路","保山道","复康路","广开二马路","广开四马路","红旗路","黄河道","密云路","密云一支路","南丰路","南京路","南开二纬路","南开三马路","三潭路","双峰道","卫津路","西湖道","西市大街","咸阳路","雅安道","延安路","玉泉路")#roads=c("广开二马路","广开四马路","红旗路","黄河道","密云路","密云一支路","南丰路","南京路","南开二纬路","南开三马路","三潭路","双峰道","卫津路","西湖道","西市大街","咸阳路","雅安道","延安路","玉泉路")#files=paste("E:/R-WorkPath/data/04",filenum,sep="","/",road,".csv")H=c();#所有的H存放地for(road in roads){files=paste("I:/data/04",filenum,sep="","/",road,".csv")GPSH=c();for(file in files){table=read.csv(file)GPS=table$gps #GPS：GPS速度（速度list）E=c() #期望E(R/S)for(i in 1:96){if(!is.numeric(GPS[i])&&!plex(GPS[i])&&is.logical(GPS[i])){waring("argument is not numeric or logical：returning NA");return(NA_real_)}}if(!is.numeric(GPS[96])&&!plex(GPS[96])&&!is.logical(GPS[96])){co=c(0,0,0,0,0,0,0,0,0,0,0);GPSH=c(GPSH,co);break;}#第一步 求均值m=mean(GPS) #m:GPS速度均值（num）#第二步 求平均调整Y=GPS-m #GPS速度偏差（list）#第三步 求累计偏离Z=c();for(i in 1:96){ Z[i]=sum(Y[1:i])}#第四步 计算值域R =c(); #求Range(1到某个数) 96个 for(i in 1:96){ R[i]=max(Z[1:i])-min(Z[1:i]) }#第五步 求标准差S =c(); #求标准差1到某个数) 96个 (被除数)for(i in 1:96){ S[i]=sd(GPS[1:i]) if(is.na(S[i]))S[i]=1;if(S[i]==0)S[i]=1;}#第六步求Rn/Sn#初始化divide = R/S #求均值 96个 理论上for(i in 1:96){if(is.na(divide[i])){divide[i]=0;}else{E[i]=mean(divide[1:i]);}}#拟合#准备beg = 1;for(i in 1:96){if(E[i]==0)beg=beg+1 #如果这个元素为0的话就不要这个数字if(is.na(E[i])) beg = beg +1;}if(beg>96){GPSH=c(GPSH,0);break;}E=E[beg:96]N=c(beg:96)zz=lm(log(E)~log(N))GPSH=c(GPSH,zz$coefficients[[2]]);#png(file=paste(road,"04",filenum,sep="",".png"), bg="transparent")plot(log(N),log(E))lines(log(N),fitted(zz))#dev.off();#dim(H)=c(11,length(filepaths))}H=c(H,GPSH)}dim(H)=c(length(files),length(roads));colnames(H)=paste(" ",roads," ");rownames(H)=paste("04月",filenum,"日",sep="")opar=par(no.readonly=TRUE)par(mfrow=c(4,6))N=c(1:11)for(i in 1:23){plot(N,H[,1],ylab="hurst")title(paste(roads[i],"11天内hurst指数",sep=""))abline(h=0)}par(opar)
filenum=c(11:15,19:24)roads=c("鞍山西道","白堤路","保山道","复康路","广开二马路","广开四马路","红旗路","黄河道","密云路","密云一支路","南丰路","南京路","南开二纬路","南开三马路","三潭路","双峰道","卫津路","西湖道","西市大街","咸阳路","雅安道","延安路","玉泉路")#roads=c("广开二马路","广开四马路","红旗路","黄河道","密云路","密云一支路","南丰路","南京路","南开二纬路","南开三马路","三潭路","双峰道","卫津路","西湖道","西市大街","咸阳路","雅安道","延安路","玉泉路")#files=paste("E:/R-WorkPath/data/04",filenum,sep="","/",road,".csv")QHMatrix_H=c();#所有的H存放地for(road in roads){files=paste("I:/data/04",filenum,sep="","/",road,".csv")QHH=c();for(file in files){table=read.csv(file)QH=table$gps #QH：QH速度（速度list）E=c() #期望E(R/S)for(i in 1:96){if(!is.numeric(QH[i])&&!plex(QH[i])&&is.logical(QH[i])){waring("argument is not numeric or logical：returning NA");return(NA_real_)}}if(!is.numeric(QH[96])&&!plex(QH[96])&&!is.logical(QH[96])){co=c(0,0,0,0,0,0,0,0,0,0,0);QHH=c(QHH,co);break;}#第一步 求均值m=mean(QH) #m:QH速度均值（num）#第二步 求平均调整Y=QH-m #QH速度偏差（list）#第三步 求累计偏离Z=c();for(i in 1:96){ Z[i]=sum(Y[1:i])}#第四步 计算值域R =c(); #求Range(1到某个数) 96个 for(i in 1:96){ R[i]=max(Z[1:i])-min(Z[1:i]) }#第五步 求标准差S =c(); #求标准差1到某个数) 96个 (被除数)for(i in 1:96){ S[i]=sd(QH[1:i]) if(is.na(S[i]))S[i]=1;if(S[i]==0)S[i]=1;}#第六步求Rn/Sn#初始化divide = R/S #求均值 96个 理论上for(i in 1:96){if(is.na(divide[i])){divide[i]=0;}else{E[i]=mean(divide[1:i]);}}#拟合#准备beg = 1;for(i in 1:96){if(E[i]==0)beg=beg+1 #如果这个元素为0的话就不要这个数字if(is.na(E[i])) beg = beg +1;}if(beg>96){QHH=c(QHH,0);break;}E=E[beg:96]N=c(beg:96)zz=lm(log(E)~log(N))QHH=c(QHH,zz$coefficients[[2]]);#png(file=paste(road,"04",filenum,sep="",".png"), bg="transparent")plot(log(N),log(E))lines(log(N),fitted(zz))#dev.off();#dim(QHMatrix_H)=c(11,length(filepaths))}QHMatrix_H=c(QHMatrix_H,QHH)}dim(QHMatrix_H)=c(length(files),length(roads));colnames(QHMatrix_H)=paste(" ",roads," ");rownames(QHMatrix_H)=paste("04月",filenum,"日",sep="")opar=par(no.readonly=TRUE)par(mfrow=c(4,6))N=c(1:11)for(i in 1:23){plot(N,QHMatrix_H[,1],ylab="hurst")title(paste(roads[i],"11天内hurst指数",sep=""))abline(h=0)}par(opar)
opar=par(no.readonly=TRUE)par(mfrow=c(4,6))N=c(1:11)a=c()for(i in 1:length(roads)){plot(H[,i],QHMatrix_H[,i],xlab="GPS Hurst 指数",ylab="切换 Hurst 指数")pz=lm(H[,i]~QHMatrix_H[,i])lines(H[,i],fitted(pz),col="red",pch=3)a=c(a,pz$coefficients[[1]])}par(opar)
