[R] Shading area under PDF of t -distribution

[Dieter Menne]

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,
            ...))
    }
}