Hey guys. I got bored and a thread reminded me that I wanted to make a function to do generic shading of distribution functions. The code could probably be a little bit better but hopefully this could help a few people out.
Code:shadeplot <- function(func = dnorm, xlim = c(-3,3), ylim = NULL, shadelim = c(2,Inf), col = "red", xlab = "", ylab = "", main = "", ...) { ## func : The density function to be plotted ## xlim : The x-limits for the plotting region ## ylim : The y-limits for the plotting region ## - If left as NULL then the default height ## - will be used and the lower limit will be ## - set to 0 ## - Testing has not been on anything other than ## - ylim = c(0, some_other_maximum_value) ## shadlim : The area that you want shaded ## col : The color to shade the region ## xlab : Label for the x-axis ## ylab : Label for the y-axis ## main : Plot title ## ... : Additional parameters to be passed on to func xs <- seq(xlim[1], xlim[2], length.out = 1000) ys <- func(xs, ...) # Find appropriate shading region lb <- max(shadelim[1], xlim[1]) ub <- min(shadelim[2], xlim[2]) # Make sure plot has lower bound of y = 0 xy <- xy.coords(xs, ys) if (is.null(ylim)){ ylim <- range(xy$y[is.finite(xy$y)]) ylim[1] <- 0 } # Plot the function plot(xs, ys, type = "l", xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, main = main) # Find which xs are inside of our bounds idx <- ((xs >= lb) & (xs <= ub)) # Shade the appropriate region. polygon(c(lb, lb, xs[idx], ub, ub), c(0, func(lb, ...), ys[idx], func(ub, ...), 0), col = col) }
Edit: On second thought it would probably be better to just create a function to shade a plot after the fact. Similar to adding points or lines... Maybe I'll do that later. That way I could remove some of the unnecessary plotting parameters (xlab, ylab, main) and it would follow a more unix-y philosophy (which is great).Code:################################# ### Testing: Examples follow ################################# # Defaults shadeplot() # F distribution shadeplot(df, xlim = c(0,5), shadelim = c(2.4, Inf), df1 = 4, df2 = 20) # Watch a t converge to a normal for(i in 1:50){ shadeplot(dt, xlim = c(-3,3), shadelim = c(1.96,Inf), df = i, main = paste("DF: ",i), ylim = c(0, .45), xlab = paste("P(T > 2) = ", round(1-pt(2, df = i), 4)) ) Sys.sleep(.08) }
Second Edit: I renamed the thread because technically this can shade things other than pdfs. However, it will always set the lower y-limit to 0 so that could be an issue. So as long as you're working with a non-negative function (or are at least only plotting a non-negative section) it should be fine.
Last edited by Dason; 02-24-2011 at 03:57 PM.
BUMP. Mainly because I didn't have this code saved in my dropbox for some reason and I just searched the entire "share your code" thread looking for this guy.
I don't have emotions and sometimes that makes me very sad.
Didn't see this. This could be a useful teaching tool. Dason could you post the link in the share your functions and code so the record is easier to find?
"If you torture the data long enough it will eventually confess."
-Ronald Harry Coase -
I do get an error here: > shadeplot(df, xlim = c(0,5), shadelim = c(2.4, Inf), df1 = 4, df2 = 20)
Code:> for(i in 1:50){ + shadeplot(dt, xlim = c(-3,3), + shadelim = c(1.96,Inf), df = i, + main = paste("DF: ",i), ylim = c(0, .45), + xlab = paste("P(T > 2) = ", round(1-pt(2, df = i), 4)) + ) + Sys.sleep(.08) + } Error in shadeplot(dt, xlim = c(-3, 3), shadelim = c(1.96, Inf), df = i, : could not find function "func"
"If you torture the data long enough it will eventually confess."
-Ronald Harry Coase -
If I remember to tomorrow then I will. But what I really should do is add it to Dmisc
I don't have emotions and sometimes that makes me very sad.
No. Not sure what's up. I'll try running a vanilla R. This fixes it though:
where I add the func <- match.fun for your codeCode:func <- match.fun(func) ys <- func(xs, ...)
"If you torture the data long enough it will eventually confess."
-Ronald Harry Coase -
Must be something I have assigned to something in my global envir. In a vanilla session it runs fine.
"If you torture the data long enough it will eventually confess."
-Ronald Harry Coase -
One could also use clip function for this.
Dason, Thanks for your idea. I am using Your example in my presentation.Code:for(i in 1:50){ xxlim <- seq(-3,3,length.out=1000) plot(xxlim,dt(x=xxlim,df=i),xlim = c(-3,3),type="l", main = paste("DF: ",i), ylim = c(0, .45), xlab = paste("P(T > 2) = ",round(1-pt(2, df = i), 4))) clip(x1=2,x2=3,y1=0,y2=0.45) polygon(xxlim,dt(x=xxlim,df=i),xlim = c(-3,3),type="l", main = paste("DF: ",i), ylim = c(0, .45), xlab = paste("P(T > 2) = ",round(1-pt(2, df = i), 4)), col="red",add=TRUE ) Sys.sleep(.08) }
In the long run, we're all dead.
Tweet |