# How to display multiple density or distribution functions on a single plot?

What would be the easiest way to plot in the same graph, the probability density function or the cumulative distribution function of a distribution, for 3-4 different parameter values.

Let's take for instance the Weibull distribution W(a,b). I would want to plot the PDF for a=1/2,b=20 , a=1/3,b=20, a=1/4,b=20, a=1/5,b=20 in the same graph and in different colors.

What is the easiest way to do this? You can recommend any kind of software. Also, the program should let me define custom distributions.

-

The question asks for "easiest." Interpreting that in terms of either (i) lines of code, (ii) naturality of expression, or (iii) raw capabilities, I find the Mathematica solutions to be well worth considering.

For example,

Plot[Evaluate[
PDF[WeibullDistribution[#, 20]][x] & /@ {1/2, 1/3, 1/4, 1/5}], {x, 0, 1},
AxesOrigin -> {0, 0}]


produces the example in the question

and

gMixture[x_, weights_, shapes_, scales_] :=
MapThread[PDF[GammaDistribution[##]][x] &, {shapes, scales}] . weights / Total[weights];
Plot[gMixture[x, {1, 2, 3}, {2, 3, 10}, {1, 1, 1}], {x, 0, 20}, AxesOrigin -> {0, 0}]


shows what it takes to define and plot a new distribution (here, a mixture of gammas):

Need something more exotic? It's likely already part of Mathematica. E.g., here is a PDF obtained from a Jacobi theta function by normalizing its area to unity:

With[{c = NIntegrate[EllipticTheta[1, z, 1/2], {z, 0, Pi}]},
Plot[EllipticTheta[1, z, 1/2] / c, {z, 0, Pi}, Filling -> Axis]]


-
 This is very nice! Could you point me to a good documentation for the "Evaluate" function? The Mathematica "Help" doesn't tell me alot about all the options this function offers. What if I would want to use a vector of values for the second (instead of 20), as well as for the first parameter in the Weibull distribution? How do I do this? More specifically, how would I plot for instance the PDF of WeibullDistribution[1/5, 50], WeibullDistribution[1/4, 40], WeibullDistribution[1/3, 30], WeibullDistribution[1/2, 20] in the same graph? – Chris May 1 '12 at 0:36 (1) For the use of Evaluate, see the examples for Plot. (2) Emulate this: Plot[Evaluate[PDF[WeibullDistribution @@ #][t] & /@ {{1/5, 50}, {1/4, 40}, {1/3, 30}}], {t, 0, 1}] (3) For tougher questions, please visit our Mathematica site. – whuber♦ May 1 '12 at 6:28

I love R, easy and free. Here's an example:

# The par removes the "padding" from the axis
par(xaxs="i", yaxs="i")

# Initiate the x, a small "by" is neat for a smooth curve
# Can't use 0 since it gives produces an integrate() error
x <- seq(0.0001, 3, by=.01)

# Just some vanity - adding a little color :-), heat.colors(5) could be an option
colors <- c("darkred", "red", "orange", "gold", "yellow")

plot(x, type="n", ylim=c(0,3), ylab="Density")
for(i in 1:5){
lines(x, dweibull(x, shape=1/i), col=colors[i])
}
title("Weibull tests")


Gives this:

### Update

I've played around with Peter Flom's suggestion with the integrate function. The prob. function, same as above:

plot(x, type="n", ylim=c(0,1), xlim=range(x), ylab="Prob")
for(i in 1:5){
lines(x, pweibull(x, shape=1/i), col=colors[i])
}
title("Using the pweibull funciton")


Give this graph:

When using the integrate function to get the "same" graph the code looks like this:

plot(x, type="n", ylim=c(0,1), xlim=c(0, max(x)), ylab="Density")
for(i in 1:5){
t <- apply(matrix(x), MARGIN=1, FUN=function(x)
integrate(function(a) dweibull(a, shape=1/i), 0, x)\$value)

lines(x, t, col=colors[i])
}
title("Using the integrate funciton")


And this gives virtually an identical graph:

-
Hi! That looks very nice, thank you very much. Is there a quick way to insert a custom function as a density function in your R code? – Chris Apr 29 '12 at 21:10
If you can write it as a function, you can probably write it in R. If (as seems likely) the pdf involves integration, have a look at this, or try ?integrate from R. – Peter Flom Apr 29 '12 at 22:04
@PeterFlom: Thank you for your tip. Added an example based on this. – Max Gordon Apr 30 '12 at 7:32

I like R too. Here is a more or less generic function to plot any probability distribution from the base R functions. It should not be difficult to extend the code with functions available in other packages, e.g. SuppDists.

plot.func <- function(distr=c("beta", "binom", "cauchy", "chisq",
"exp", "f","gamma", "geom", "hyper",
"logis", "lnorm", "nbinom", "norm",
"pois", "t", "unif", "weibull"),
what=c("pdf", "cdf"), params=list(), type="b",
xlim=c(0, 1), log=FALSE, n=101, add=FALSE, ...) {
what <- match.arg(what)
d <- match.fun(paste(switch(what, pdf = "d", cdf = "p"),
distr, sep=""))
# Define x-values (because we won't use 'curve') as last parameter
# (with pdf, it should be 'x', while for cdf it is 'q').
len <- length(params)
params[[len+1]] <- seq(xlim[1], xlim[2], length=n)
if (add) lines(params[[len+1]], do.call(d, params), type, ...)
else plot(params[[len+1]], do.call(d, params), type, ...)
}


It's a bit crappy and I haven't tested it a lot. The params list must obey R's conventions for naming {C|P}DF parameters (e.g., shape and scale for the Weibull distribution, and not a or b). There's room for improvement, especially about the way it handles multiple plotting on the same graphic device (and, actually, passing vector of parameters only works as a side-effect when type="p"). Also, there's not much parameter checking!

Here are some examples of use:

# Normal CDF
xl <- c(-5, 5)
plot.func("norm", what="pdf",  params=list(mean=1, sd=1.2),
xlim=xl, ylim=c(0,.5), cex=.8, type="l", xlab="x", ylab="F(x)")
plot.func("norm", what="pdf",  params=list(mean=3, sd=.8),
xlim=xl, add=TRUE, pch=19, cex=.8)
plot.func("norm", what="pdf",  params=list(mean=.5, sd=1.3), n=201,
xlim=xl, add=TRUE, pch=19, cex=.4, type="p", col="steelblue")
title(main="Some gaussian PDFs")

# Standard normal PDF
plot.func("norm", "cdf", xlab="Quantile (x)", ylab="P(X<x)", xlim=c(-3,3), type="l",
main="Some gaussian CDFs")
plot.func("norm", "cdf", list(sd=c(0.5,1.5)), xlim=c(-3,3), add=TRUE,
type="p", pch=c("o","+"), n=201, cex=.8)
legend("topleft", paste("N(0;", c(1,0.5,1.5), ")", sep=""),
lty=c(1,NA,NA), pch=c(NA,"o","+"), bty="n")

# Weibull distribution
s <- c(.5,.75,1)
plot.func("weibull", what="pdf", xlim=c(0,1), params=list(shape=s),
col=1:3, type="p", n=301, pch=19, cex=.6, xlab="", ylab="")
title(main="Weibull distribution", xlab="x", ylab="F(x)")
legend("topright", legend=as.character(s), title="Shape", col=1:3, pch=19)


-
Wow, impressive code - really interesting to see how to take my idea to the next level. Thanks @chl! – Max Gordon Apr 30 '12 at 15:13
I can't quite figure out why looking over the code at first glance, but the pdf's and the cdf's are switched. – Andy W Apr 30 '12 at 20:04
@Andy Sure! It looks like I swapped the type arguments at some point when cleaning up an earlier version of the code. (And good catch, btw) – chl Apr 30 '12 at 20:26
Thank you so much! – Chris May 1 '12 at 0:32