So zeichnen Sie ein Quantilband (in R) auf
https://stackoverflow.com//questions/12658555
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianFrage
Ich habe eine CSV-Datei, die Zeilen für jedes (Java GC-)Ereignis enthält, an dem ich interessiert bin.Das Objekt besteht aus einem Zeitstempel im Sekundenbruchteil (nicht äquidistant) und einigen Variablen.Das Objekt sieht so aus:
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)
Ergebnisse in
'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 ...
In diesem Fall ist mir die Pausenzeit (in Sekunden) wichtig.Ich möchte ein Diagramm zeichnen, das mir für jede (Wanduhr-)Stunde grundsätzlich den Mittelwert als Linie, die 2 % und 98 % als grauen Korridor und den Maximalwert (innerhalb jeder Stunde) als rote Linie anzeigt.
Ich habe einige Arbeit geleistet, aber die Verwendung der q98-Funktionen ist hässlich, die Verwendung mehrzeiliger Anweisungen scheint verschwenderisch zu sein, und ich weiß nicht, wie ich eine Grauzone zwischen q02 und q98 erreichen kann:
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")
Dies ergibt nun ein Diagramm, das schwarze Punkte als Maximum, eine blaue Linie als Stundendurchschnitt und eine untere und obere grüne Linie von 0,2 + 0,98 aufweist.Ich denke, es wäre besser lesbar, einen grauen Korridor zu haben, vielleicht eine gestrichelte Maximumlinie (rot) und die Achsenbeschriftungen irgendwie zu korrigieren.
Irgendwelche Vorschläge?(die Datei ist oben verfügbar)
Lösung
Schön, hier alte Debian-Kollegen zu sehen :) Deine Antwort ist schon ziemlich nett.Da ich zufällig viel mit Zeitreihen arbeite, dachte ich, ich füge eine Variante hinzu, die das Excellent verwendet Zoo Und xts Pakete.Letzteres baut auf ersterem auf und verfügt unter anderem über die period.apply() Funktion, die wir hier zusammen mit der verwenden können endpoints() Funktion, mit der Sie zweistündliche Aggregate erhalten.
Also oben würde ich verwenden
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
eine erstellen zoo Objekt.Ihre Funktion bleibt bestehen:
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)
}
Und dann können wir es etwas vereinfachen:
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()
was gibt
Die Achse ist jetzt automatisch und zeigt nur Wochentage an, da die Spanne weniger als eine Woche beträgt;Dies kann bei Bedarf überschrieben werden.
Andere Tipps
Das musst du ausprobieren polygon.Dieser Code kann nützlich sein:
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)
oder
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)
wo die Funktion plotAreaCI ist definiert durch:
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())
}
Dies ist der Code, den ich verwende, um die zeitliche Variation von Laboranalyten (in diesem Fall des systolischen Blutdrucks) darzustellen:
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")
Sollte nicht allzu schwer sein, das an Ihr Problem anzupassen ...oder Sie könnten ein reproduzierbares Beispiel posten, mit dem Sie arbeiten können.Dies ergibt das 10., 25., 50., 75., 90., 95. und 97,5. Perzentil in einem einzigen Datenrahmen und Matplot übernimmt die Darstellung eines solchen Objekts.
Die Grauzone?, ...Der übliche Ansatz besteht darin, ein Polygon zu zeichnen, das an den unteren Grenzen ausgeht, sich am rechten Rand „dreht“, auf der oberen Seite zurückkommt und auf der linken Seite wieder eine Verbindung herstellt.Der polygon Argumente werden als eingerichtet x, y.Da ist ein col Argument, das Sie auf „grau“ setzen würden.
Um „2-Stunden“-Sequenzen zu erstellen, mit denen Sie Ihren Datenrahmen zusammenführen oder verwenden könnencut.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"
Ich habe es nicht dokumentiert gesehen, aber Sie können Vielfache von Intervallen verwenden 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 ...
Ich bin derzeit nicht im folgenden Skript angekommen (muss immer noch die fortgeschrittenere Antwort von Dwin ansehen).Es sieht jetzt etwas aus, als hätte ich gesucht, aber der Code ist immer noch ziemlich hässlich (zum Beispiel weiß ich nicht, wie ich das Etikett ausrichten kann und wie man ordnungsgemäße XLAB-Etiketten ermittelt): generasacodicetagpre.
