クォンティル バンドをプロットする方法 (R 単位)
https://stackoverflow.com//questions/12658555
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian質問
興味のある各 (Java GC) イベントの行を含む CSV ファイルがあります。オブジェクトは、1 秒未満のタイムスタンプ (非等距離) といくつかの変数で構成されます。オブジェクトは次のようになります。
gcdata <- read.table("http://bernd.eckenfels.net/view/gc1001.ygc.csv",header=TRUE,sep=",", dec=".")
start = as.POSIXct(strptime("2012-01-01 00:00:00", format="%Y-%m-%d %H:%M:%S"))
gcdata.date = gcdata$Timestamp + start
gcdata = gcdata[,2:7] # remove old date col
gcdata=data.frame(date=gcdata.date,gcdata)
str(gcdata)
での結果
'data.frame': 2997 obs. of 7 variables:
$ date : POSIXct, format: "2012-01-01 00:00:06" "2012-01-01 00:00:06" "2012-01-01 00:00:18" ...
$ Distance.s. : num 0 0.165 11.289 9.029 11.161 ...
$ YGUsedBefore.K.: int 1610619 20140726 20148325 20213304 20310849 20404772 20561918 21115577 21479211 21544930 ...
$ YGUsedAfter.K. : int 7990 15589 80568 178113 272036 429182 982841 1346475 1412181 1355412 ...
$ Promoted.K. : int 0 0 0 0 8226 937 65429 71166 62548 143638 ...
$ YGCapacity.K. : int 22649280 22649280 22649280 22649280 22649280 22649280 22649280 22649280 22649280 22649280 ...
$ Pause.s. : num 0.0379 0.022 0.0287 0.0509 0.109 ...
この場合、一時停止時間 (秒単位) が重要です。図をプロットしたいと思います。これは、（壁時計の）時間ごとに、基本的に平均を線で、2% と 98% を灰色の回廊で、(各時間内の) 最大値を赤い線で表示します。
いくつかの作業を行いましたが、q98 関数の使用は醜く、複数行のステートメントを使用しなければならないのは無駄なようで、q02 と q98 の間の灰色の領域を実現する方法がわかりません。
q02 <- function(x, ...) { x <- quantile(x,probs=c(0.2)) }
q98 <- function(x, ...) { x <- quantile(x,probs=c(0.98)) }
hours = droplevels(cut(gcdata$date, breaks="hours")) # can I have 2 hours?
plot(aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=max),ylim=c(0,2), col="red", ylab="Pause(s)", xlab="Days") # Is always black?
lines(aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=q98),ylim=c(0,2), col="green")
lines(aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=q02),ylim=c(0,2), col="green")
lines(aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=mean),ylim=c(0,2), col="blue")
これにより、黒い点が最大値、青い線が時間平均、下限と上限が 0,2 + 0,98 の緑色の線を持つチャートが作成されます。灰色のコリドー、おそらく最大値の破線 (赤色) を表示し、何らかの方法で軸ラベルを修正した方が読みやすくなると思います。
助言がありますか？(ファイルは上記から入手可能です)
解決
ここで仲間のDebianの古いタイマーを見るのが良い:)あなたの答えはすでにかなり素敵です。私が時系列でたくさん作業するようになるので、私は優れたperiod.apply()関数が2時間の集計を得るためにここで使用できる。
だから私はを使いたいです
library(zoo) # for zoo objects
library(xts) # for period.apply
gcdata <- read.table("http://bernd.eckenfels.net/view/gc1001.ygc.csv",
header=TRUE, sep=",", dec=".")
timestamps <- gcdata$Timestamp +
as.POSIXct(strptime("2012-01-01 00:00:00",
format="%Y-%m-%d %H:%M:%S"))
gcdatazoo <- zoo(gcdata[-1], order.by=timestamps) # as zoo object
endpoints()オブジェクトを作成します。あなたの機能は残ります：
plotAreaCorridor <- function(x, y, col.poly1="lightgray", col.poly2="gray",...) {
x.pol <- c(x, rev(x), x[1])
y.pol <- c(y[,1], rev(y[,5]),y[,1][1])
plot(x, y[,6]+1, type="n", ...)
polygon(x.pol, y.pol, col=col.poly1, lty=0)
x.pol <- c(x, rev(x), x[1])
y.pol <- c(y[,2], rev(y[,4]), y[,1][1])
polygon(x.pol, y.pol, col=col.poly2, lty=0)
lines(x, y[,3], col="blue") # median
lines(x, y[,6], col="red") # max
invisible(NULL)
}
そして私達はそれから少し単純化することができます：
agg <- period.apply(gcdatazoo[,"Pause.s."], # to which data
INDEX=endpoints(gcdatazoo, "hours", k=2), # every 2 hours
FUN=function(x) quantile(x, # what fun.
probs=c(5,20,50,80,95,100)/100))
#v99 = q99(gcdata$Pause.s.) # what is q99 ?
v99 <- mean(agg[,5]) # mean of 95-th percentile?
plotAreaCorridor(index(agg), # use time index as x axis
coredata(agg), # and matrix part of zoo object as data
ylim=c(0,max(agg[,5])*1.5),
ylab="Quantiles of GC events",
main="NewPar Collection Activity")
abline(h=median(gcdatazoo[,"Pause.s."]), col="lightblue")
abline(h=v99, col="grey")
labeltxt <- paste("99%=",round(v99,digits=3),"s n=", nrow(gcdatazoo),sep="")
text(x=index(agg)[20], y=1.5*v99, labeltxt, col="grey", pos=3) # or legend()
軸は自動的に、スパンが週未満の場合にのみ平日を示します。これは必要に応じて上書きすることができます。
他のヒント
polygonを試してみてください。このコードは便利です：
y98 = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=q98)
y02 = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=q02)
ymax = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=max)
ymin = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=min)
ymean = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=mean)
x = ymean[,1]
y1 = cbind(y02[,2], ymean[,2], y98[,2])
y2 = cbind(ymin[,2], ymean[,2], ymax[,2])
plotAreaCI(x,y2, ylim=c(0,2), xlab="time", ylab="variable")
plotAreaCI(x,y1, ylim=c(0,2), poly.col="blue", add=TRUE)
または
plotAreaCI(x,y2, ylim=c(0,2), xlab="time", ylab="variable", nice.x = TRUE)
plotAreaCI(x,y1, ylim=c(0,2), mean.lwd=2, poly.col="blue", add=TRUE)
関数plotAreaCIが次のように定義されている場合：
plotAreaCI = function(x, y, add=FALSE, nice.x = FALSE,
xlim=NULL, ylim=NULL,
mean.col="black", mean.lwd=1.5,
poly.col="gray", poly.lty=3,
xlab=NULL, ylab=NULL, main="",
...) {
isFactorX = isClass("factor", x)
if(isFactorX) {
x.label = x
x = as.numeric(x)
}
if(is.null(xlim)) xlim=range(x, na.rm=TRUE)
if(is.null(ylim)) ylim=range(y, na.rm=TRUE)
x.pol = c(x, rev(x), x[1])
y.pol = c(y[,1], rev(y[,3]), y[,1][3])
if(!add) {
plot.new()
plot.window(xlim=xlim, ylim=ylim, ...)
if(!nice.x & isFactorX) {
axis(1, at=x, labels=x.label)
} else {
xticks = axTicks(1)
if(isFactorX) {
xticks = xticks[xticks>=1]
axis(1, at=xticks, labels=x.label[xticks])
} else {
axis(1)
}
}
axis(2, las=1)
box()
title(xlab=xlab, ylab=ylab, main=main)
}
polygon(x.pol, y.pol, col=poly.col, lty=poly.lty)
lines(x, y[,2], col=mean.col, lwd=mean.lwd)
return(invisible())
}
これは、実験室分析対象物 (この場合は最高血圧) の時間的変化をプロットするために使用するコードです。
SBP.qtr.mat <- aggregate(set1HLI$SBP,
list( year(set1HLI$Drawdt)+0.25* quarter(set1HLI$Drawdt)),
quantile, prob=c(0.1,0.25,0.5,0.75, 0.9,0.95, 0.975), na.rm=TRUE)
matplot(SBP.qtr.mat[,1], SBP.qtr.mat$x, type="pl")
それをあなたの問題に適応させるのはそれほど難しいことではないはずです...または、使用できる再現可能な例を投稿することもできます。これにより、単一の data.frame 内の 10、25、50、75、90、95、および 97.5 パーセンタイルが得られ、matplot はそのようなオブジェクトのプロットを処理します。
グレーゾーン?、...通常のアプローチは、下限から出て右端で「回転」し、高い側に戻り、左側で再び接続する多角形をプロットすることです。の polygon 引数は次のように設定されます x, y. 。があります col 引数を「グレー」に設定します。
データフレームをマージしたり、cut.POSIXt" as a breaks argument , there is the option of using multiples of time units withseq.POSIXt`:
> seq(ISOdate(1910,1,1), ISOdate(1999,1,1), "10 years")
[1] "1910-01-01 12:00:00 GMT" "1920-01-01 12:00:00 GMT" "1930-01-01 12:00:00 GMT" "1940-01-01 12:00:00 GMT"
[5] "1950-01-01 12:00:00 GMT" "1960-01-01 12:00:00 GMT" "1970-01-01 12:00:00 GMT" "1980-01-01 12:00:00 GMT"
[9] "1990-01-01 12:00:00 GMT"
文書化されていないのですが、interval の倍数を使用できます。 cut.POSIXt:
> str( cut( seq(ISOdate(1910,1,1), ISOdate(1999,1,1), "years"), "10 years") )
Factor w/ 9 levels "1910-01-01","1920-01-01",..: 1 1 1 1 1 1 1 1 1 1 ...
> str( cut( seq(ISOdate(1910,1,1), ISOdate(1999,1,1), "years"), "5 years") )
Factor w/ 18 levels "1910-01-01","1915-01-01",..: 1 1 1 1 1 2 2 2 2 2 ...
現在次のスクリプトに到着していません（まだDWINからのより高度な回答を見る必要があります）。私が探していたようなようになりましたが、コードはまだかなり醜いです（たとえば、ラベルを揃える方法と適切なXlabラベルの取得方法を知らない）：
plotAreaCorridor = function(x, y, col.poly1="lightgray", col.poly2="gray",...) {
x.pol = c(x, rev(x), x[1])
y.pol = c(y[,1], rev(y[,5]),y[,1][1])
plot(x, y[,6]+1, type="n", ...) # ugly since type="n" does not work for factor
polygon(x.pol, y.pol, col=col.poly1, lty=0)
x.pol = c(x, rev(x), x[1])
y.pol = c(y[,2], rev(y[,4]), y[,1][1])
polygon(x.pol, y.pol, col=col.poly2, lty=0)
lines(x, y[,3], col="blue") # median
lines(x, y[,6], col="red") # max
return(invisible())
}
pause = gcdata$Pause.s.
hours = droplevels(cut(gcdata$date, breaks="hours")) # can I have 2 hours?
agg = aggregate(pause ~ hours, FUN=quantile, probs=c(5,20,50,80,95,100)/100)
x = agg$hours
ys = agg$pause
q99 <- function(x, ...) { x <- quantile(x,probs=c(0.99)) }
v99 = q99(gcdata$Pause.s.)
vmed = median(gcdata$Pause.s.)
plotAreaCorridor(x, ys,ylim=c(0,v99*1.5))
abline(h=vmed, col="lightblue")
abline(h=v99, col="grey")
label=paste("99%=",round(v99,digits=3),"s n=", length(gcdata$date),sep="")
text(x=30, y=v99, label, col="grey", pos=3)
title("NewPar Collection Activity")
所属していません StackOverflow
