[R] Shading area under PDF of t -distribution
Dieter Menne
dieter.menne at menne-biomed.de
Fri Nov 2 10:26:36 CET 2007
Chung-hong Chan <chainsawtiney <at> gmail.com> writes:
> I have plot the PDF of t distribution with df = 74.
> curve(dt(x,df=74),from=-4, to=4)
>
> how can I shade the area under curve (for example, col="red") from t=+- 1.996?
If you could used lattice plots instead, below is a modified version of the
standard lattice densityplot. You can use it with a call like:
densityplot(form, data=DHL, shadelimit=shadelimit,
panel=panel.shadeddensityplot,
shadecol='#FFEECC')
I am sure Deepayan will have a more elegant solution, but this worked for me.
Dieter
------
# Densityplot with limits
panel.shadeddensityplot =
function (x, darg = list(n = 100), plot.points = "jitter", ref = FALSE,
groups = NULL, jitter.amount = 0.01 * diff(current.panel.limits()$ylim),
type = "p",shadelimit=0, shadecol="red", ...)
{
x=na.omit(x)
if (ref) {
reference.line <- trellis.par.get("reference.line")
panel.abline(h = 0, col = reference.line$col, lty = reference.line$lty,
lwd = reference.line$lwd)
}
plot.line <- trellis.par.get("plot.line")
superpose.line <- trellis.par.get("superpose.line")
if (!is.null(groups)) {
panel.superpose(x, darg = darg, plot.points = plot.points,
ref = FALSE, groups = groups, panel.groups = panel.densityplot,
jitter.amount = jitter.amount, type = type, ...)
}
else {
if (sum(!is.na(x)) > 1) {
h <- do.call("density", c(list(x = x), darg))
lim <- current.panel.limits()$xlim
id <- h$x > min(lim) & h$x < max(lim)
panel.lines(x = h$x[id], y = h$y[id], ...)
idp = h$x> min(lim) & h$x < shadelimit
xp = h$x[idp]
xp = c(xp[1],xp,xp[length(xp)])
yp = h$y[idp]
yp = c(0,yp,0)
panel.polygon(x = xp, y = yp,col=shadecol, ...)
grid.text(x=0.01,y=0.99,default.units="npc",
gp=gpar(fontsize=10),label=
paste("n = ",length(x[x<=shadelimit]),"/",length(x),sep=""),
just=c("left","top"))
}
switch(as.character(plot.points), "TRUE" = panel.xyplot(x = x,
y = rep(0, length(x)), type = type, ...), rug = panel.rug(x = x,
start = 0, end = 0, x.units = c("npc", "native"),
type = type, ...), jitter = panel.xyplot(x = x, y = jitter(rep(0,
length(x)), amount = jitter.amount), type = type,
...))
}
}
More information about the R-help
mailing list