¿Cómo mostrar varias funciones de densidad o distribución en un solo gráfico?

Question

Cuál sería la forma más sencilla de representar en un mismo gráfico la función de densidad de probabilidad o la función de distribución acumulativa de una distribución, para 3-4 valores diferentes de los parámetros.

Tomemos por ejemplo la distribución Weibull W(a,b). Me gustaría trazar la PDF para a=1/2,b=20 , a=1/3,b=20, a=1/4,b=20, a=1/5,b=20 en el mismo gráfico y en diferentes colores.

¿Cuál es la forma más fácil de hacerlo? Puede recomendar cualquier tipo de software. Además, el programa debería permitirme definir distribuciones personalizadas.

Answer 1

Me encanta R fácil y gratuito. He aquí un ejemplo:

# 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")

Da esto:

Actualización

He jugado con la sugerencia de Peter Flom con la función integrar. La función prob., igual que la anterior:

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

Da esta gráfica:

Cuando se utiliza la función de integración para obtener la "misma" gráfica el código se ve así:

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

Y esto da un gráfico prácticamente idéntico:

Answer 2

A mí también me gusta R. Aquí hay una función más o menos genérica para graficar cualquier distribución de probabilidad del funciones base de R . No debería ser difícil ampliar el código con funciones disponibles en otros paquetes, por ejemplo 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, ...)
}

Es un poco cutre y no lo he probado mucho. El params debe obedecer las convenciones de R para nombrar los parámetros {C|P}DF (por ejemplo, shape y scale para la distribución de Weibull, y no a o b ). Hay margen de mejora, especialmente sobre la forma en que maneja el trazado múltiple en el mismo dispositivo gráfico (y, en realidad, pasar el vector de parámetros sólo funciona como un efecto secundario cuando type="p" ). Además, ¡no hay mucha comprobación de parámetros!

Estos son algunos ejemplos de uso:

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

Answer 3

La pregunta pide "lo más fácil". Interpretando eso en términos de (i) líneas de código, (ii) naturalidad de la expresión, o (iii) capacidades brutas, encuentro que el Mathematica soluciones que vale la pena considerar.

Por ejemplo,

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

produce el ejemplo de la pregunta

y

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

muestra lo que se necesita para definir y trazar una nueva distribución (aquí, una mezcla de gammas):

¿Necesitas algo más exótico? Es probable que ya forme parte de Mathematica. Por ejemplo, aquí hay un PDF obtenido de un Función theta de Jacobi normalizando su área a la unidad:

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