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

Question

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.

Answer 1

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]]

Answer 2

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:

Answer 3

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)