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,

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

produces the example in the question

Weibull distributions


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):

Gamma mixture

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

Elliptic PDF

  • $\begingroup$ 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? $\endgroup$
    – Chris
    May 1 '12 at 0:36
  • $\begingroup$ (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. $\endgroup$
    – 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:

Density Weibull


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:

Prob. Weibull

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:

Prob Weibull using the integrate function

  • $\begingroup$ 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? $\endgroup$
    – Chris
    Apr 29 '12 at 21:10
  • 1
    $\begingroup$ 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. $\endgroup$
    – Peter Flom
    Apr 29 '12 at 22:04
  • $\begingroup$ @PeterFlom: Thank you for your tip. Added an example based on this. $\endgroup$
    – 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)

enter image description here

  • $\begingroup$ Wow, impressive code - really interesting to see how to take my idea to the next level. Thanks @chl! $\endgroup$
    – Max Gordon
    Apr 30 '12 at 15:13
  • 1
    $\begingroup$ I can't quite figure out why looking over the code at first glance, but the pdf's and the cdf's are switched. $\endgroup$
    – Andy W
    Apr 30 '12 at 20:04
  • $\begingroup$ @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) $\endgroup$
    – chl
    Apr 30 '12 at 20:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.