倍率パーセンタイルスクリプト
apply(mat_food,"median",2)
medi1_food<-NULL
mean1_food<-apply(mat_food,2,"mean")
for(i in 1:length(mat_food[,1])){
q<-(mat_food[i,]/mean1_food)[which(mean1_food>=100)]
medi1_food[i]<-median(q)
}
#r0<-(mean_end+10^-3)/(mean_start+10^-3)
mean_start_b<-NULL
mean_end_b<-NULL
for(i in 1:length(mat_food[1,])){
sst<-as.Date("2007-01-01")-as.Date("2006-10-31")
fft<-as.Date("2007-12-31")-as.Date("2006-10-31")
sst2<-as.Date("2014-01-01")-as.Date("2006-10-31")
fft2<-as.Date("2014-12-31")-as.Date("2006-10-31")
mean_start_b[i]<-median(mat_food[,i][sst:fft]/medi1_food[sst:fft])
mean_end_b[i]<-median(mat_food[,i][sst2:fft2]/medi1_food[sst2:fft2])
}
r0<-(mean_end_b+0)/(mean_start_b+0)
cc<-exp(seq(log(min(mean_start+10^-3,na.rm=T)+0.001),log(max(mean_start,na.rm=T)+0.001),length.out=20))
k<-1
m1<-NULL
s1<-NULL
q99<-NULL
q95<-NULL
q50<-NULL
q25<-NULL
q75<-NULL
q10<-NULL
q05<-NULL
q01<-NULL
q001<-NULL
for(i in 1:length(cc[-1])){
#for(i in 3:3){
v<-which(mean_start_b>=cc[i-1] & mean_start_b<=cc[i])
rr<-r0[v]
if(length(rr)>=30){
if(k==1){
plot(sort(log(rr)),length(rr):1/length(rr),type="l",col= topo.colors(length(cc))[i],log="xy")
m1[i]<-mean(log(rr))
s1[i]<-sd(log(rr))
q99[i]<-quantile(log(rr),0.99)
q95[i]<-quantile(log(rr),0.95)
q90[i]<-quantile(log(rr),0.9)
q50[i]<-quantile(log(rr),0.5)
q25[i]<-quantile(log(rr),0.25)
q75[i]<-quantile(log(rr),0.75)
q10[i]<-quantile(log(rr),0.1)
q05[i]<-quantile(log(rr),0.05)
q001[i]<-quantile(log(rr),0.01)
k=k+1
}else{
points(sort(log(rr)),length(rr):1/length(rr),type="l",col= topo.colors(length(cc))[i],xlim=c(0,1))
m1[i]<-mean(log(rr))
s1[i]<-sd(log(rr))
q95[i]<-quantile(log(rr),0.95)
q90[i]<-quantile(log(rr),0.9)
q99[i]<-quantile(log(rr),0.99)
q50[i]<-quantile(log(rr),0.5)
q25[i]<-quantile(log(rr),0.25)
q75[i]<-quantile(log(rr),0.75)
q10[i]<-quantile(log(rr),0.1)
q05[i]<-quantile(log(rr),0.05)
q001[i]<-quantile(log(rr),0.01)
k=k+1
}
}
}
#plot(cc[4:length(cc)],exp(m1),log="x")
qplot(cc[4:length(cc)],exp(m1),xlab="Mean of start",ylab="Mean of 6 year ratio",log="x",geom=c("line","point"))+theme_bw(base_size=30)
qplot(cc[4:length(cc)],exp(s1),xlab="Mean of Start",ylab="SD of 6 year ratio",log="x",geom=c("line","point"),ylim=c(0,15))+theme_bw(base_size=30)
qplot(cc[4:length(cc)],exp(q90),xlab="Mean of Start",ylab="SD of 6 year ratio",log="x",geom=c("line","point"),ylim=c(0,15))+theme_bw(base_size=30)
v<-qplot(cc[2:length(cc)],exp(q90),xlab="Mean of Start",ylab="SD of 6 year ratio",log="x",geom=c("line","point"),ylim=c(0,15))+theme_bw(base_size=30)
v<-v+geom_point(cc[4:length(cc)],exp(q75),xlab="Mean of Start",ylab="SD of 6 year ratio",log="x",ylim=c(0,15))+theme_bw(base_size=30)
v<-v+geom_point(cc[4:length(cc)],exp(q50),xlab="Mean of Start",ylab="SD of 6 year ratio",log="x",ylim=c(0,15))+theme_bw(base_size=30)
v<-v+geom_point(cc[4:length(cc)],exp(q25),xlab="Mean of Start",ylab="SD of 6 year ratio",log="x",ylim=c(0,15))+theme_bw(base_size=30)
v<-v+geom_point(cc[4:length(cc)],exp(q10),xlab="Mean of Start",ylab="SD of 6 year ratio",log="x",ylim=c(0,15))+theme_bw(base_size=30)
print(v)
plot(cc[3:length(cc)],exp(q95),type="b",ylim=c(0.01,20),log="xy",col=1)
points(cc[3:length(cc)],exp(q99),type="b",ylim=c(0.1,5),log="xy",col=2)
points(cc[3:length(cc)],exp(q90),type="b",ylim=c(0.1,5),log="xy",col=2)
points(cc[3:length(cc)],exp(q25),type="b",col=3)
points(cc[3:length(cc)],exp(q50),type="b",col=4)
points(cc[3:length(cc)],exp(q75),type="b",col=5)
points(cc[3:length(cc)],exp(q10),type="b",col=6)
points(cc[3:length(cc)],exp(q05),type="b",col=7)
points(cc[3:length(cc)],exp(q01),type="b",col=7)
points(cc[3:length(cc)],exp(q001),type="b",col=8)
abline(h=1.0,lty=2)
plot(cc[3:length(cc)],exp(q95),type="b",ylim=c(0.01,20),log="xy",cex=2,cex.axis=1.5,pch=19,col=4,xlab="Mean 2007",ylab="Percentile Ratio 2014/2007",cex.axis=1.5,tck=0.01)
points(cc[3:length(cc)],exp(q99),type="b",ylim=c(0.1,5),log="xy",col=5,cex=2,cex.axis=1.5,pch=7)
points(cc[3:length(cc)],exp(q90),type="b",ylim=c(0.1,5),log="xy",cex=2,cex.axis=1.5,pch=15,col=3)
points(cc[3:length(cc)],exp(q25),type="b",cex=2,cex.axis=1.5,pch=18,col=2)
points(cc[3:length(cc)],exp(q50),type="b",cex=2,cex.axis=1.5,pch=17,col=1)
points(cc[3:length(cc)],exp(q75),type="b",cex=2,cex.axis=1.5,pch=18,col=2)
points(cc[3:length(cc)],exp(q10),type="b",cex=2,cex.axis=1.5,pch=15,col=3)
points(cc[3:length(cc)],exp(q05),type="b",cex=2,cex.axis=1.5,pch=19,col=4)
#points(cc[3:length(cc)],exp(q01),type="b",cex=2,cex.axis=1.5,pch=6)
points(cc[3:length(cc)],exp(q001),type="b",col=5,cex=2,cex.axis=1.5,pch=7)
abline(h=1.0,lty=2)
