CODEWITHSUNDEEP
r
user297850
user297850

Reputation: 8015

How to plot the probabilistic density function of a function?

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

Upvotes: 7

Views: 4203

Answers (4)

Greg Snow
Greg Snow

Reputation: 49670

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

Upvotes: 9

JD Long
JD Long

Reputation: 60756

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

Upvotes: 7

nullglob
nullglob

Reputation: 7043

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)

Upvotes: 2

Jim Lewis
Jim Lewis

Reputation: 45125

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.

Upvotes: 1

Related Questions