R语言 Hurst指数计算

GPS

filenum=c(11:15,19:24)

roads=c("鞍山西道","白堤路","保山道","复康路","广开二马路","广开四马路","红旗路","黄河道","密云路","密云一支路","南丰路","南京路","南开二纬路","南开三马路","三潭路","双峰道","卫津路","西湖道","西市大街","咸阳路","雅安道","延安路","玉泉路")
#roads=c("广开二马路","广开四马路","红旗路","黄河道","密云路","密云一支路","南丰路","南京路","南开二纬路","南开三马路","三潭路","双峰道","卫津路","西湖道","西市大街","咸阳路","雅安道","延安路","玉泉路")

#files=paste("E:/R-WorkPath/data/201104",filenum,sep="","/",road,".csv")


H=c();#所有的H存放地
for(road in roads){

	files=paste("I:/data/201104",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])&&!is.complex(GPS[i])&&is.logical(GPS[i])){
				waring("argument is not numeric or logical：returning NA");
				return(NA_real_)
			}
		}
		if(!is.numeric(GPS[96])&&!is.complex(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,"201104",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("2011年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/201104",filenum,sep="","/",road,".csv")


QHMatrix_H=c();#所有的H存放地
for(road in roads){

	files=paste("I:/data/201104",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])&&!is.complex(QH[i])&&is.logical(QH[i])){
				waring("argument is not numeric or logical：returning NA");
				return(NA_real_)
			}
		}
		if(!is.numeric(QH[96])&&!is.complex(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,"201104",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("2011年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)