Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to plot the probabilistic density function of a function?

Tags:

r

Assume A follows Exponential distribution; B follows Gamma distribution How to plot the PDF of 0.5*(A+B)

like image 472
user297850 Avatar asked Aug 25 '10 04:08

user297850


4 Answers

This is fairly straight forward using the "distr" package:

library(distr)

A <- Exp(rate=3)
B <- Gammad(shape=2, scale=3)

conv <- 0.5*(A+B)

plot(conv)
plot(conv, to.draw.arg=1)

Edit by JD Long

Resulting plot looks like this: alt text

like image 115
Greg Snow Avatar answered Nov 03 '22 13:11

Greg Snow


If you're just looking for fast graph I usually do the quick and dirty simulation approach. I do some draws, slam a Gaussian density on the draws and plot that bad boy:

numDraws   <- 1e6
gammaDraws <- rgamma(numDraws, 2)
expDraws   <- rexp(numDraws)
combined   <- .5 * (gammaDraws + expDraws)
plot(density(combined))

output should look a little like this:

alt text

like image 25
JD Long Avatar answered Nov 03 '22 14:11

JD Long


Here is an attempt at doing the convolution (which @Jim Lewis refers to) in R. Note that there are probably much more efficient ways of doing this.

lower <- 0
upper <- 20
t <- seq(lower,upper,0.01)
fA <- dexp(t, rate = 0.4)
fB <- dgamma(t,shape = 8, rate = 2)
## C has the same distribution as (A + B)/2
dC <- function(x, lower, upper, exp.rate, gamma.rate, gamma.shape){
  integrand <- function(Y, X, exp.rate, gamma.rate, gamma.shape){
    dexp(Y, rate = exp.rate)*dgamma(2*X-Y, rate = gamma.rate, shape = gamma.shape)*2
  }
  out <- NULL
  for(ix in seq_along(x)){
    out[ix] <-
      integrate(integrand, lower = lower, upper = upper,
                X = x[ix], exp.rate = exp.rate,
                gamma.rate = gamma.rate, gamma.shape = gamma.shape)$value
  }
  return(out)
}
fC <- dC(t, lower=lower, upper=upper, exp.rate=0.4, gamma.rate=2, gamma.shape=8)
## plot the resulting distribution
plot(t,fA,
     ylim = range(fA,fB,na.rm=TRUE,finite = TRUE),
     xlab = 'x',ylab = 'f(x)',type = 'l')
lines(t,fB,lty = 2)
lines(t,fC,lty = 3)
legend('topright', c('A ~ exp(0.4)','B ~ gamma(8,2)', 'C ~ (A+B)/2'),lty = 1:3)
like image 25
nullglob Avatar answered Nov 03 '22 14:11

nullglob


I'm not an R programmer, but it might be helpful to know that for independent random variables with PDFs f1(x) and f2(x), the PDF of the sum of the two variables is given by the convolution f1 * f2 (x) of the two input PDFs.

like image 31
Jim Lewis Avatar answered Nov 03 '22 13:11

Jim Lewis